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ABSTHACr 


A modal to pradict tha distribution of water quality 
parameters in three dimensions has baen developad, and 
applied to several theoretical and real problems. The 
mass transport equation is solved using a non- 
dimensional vertical axis and an alternating-iireotion- 
implicit finite difference technique. The reaction 
kinetics of the constituents are incorporated into a 
matrix method which permits computation of tha 
interactions of multiple constituents. 

Extensive literatura reviews were made to determine 
the most appropriate methods available for tha 
computation of dispersion coefficients and colifocm 
bacteria decay rates. Numerical investigations of 
dispersive and dissipative effects showed that the 
three-dimensional model performed as predicted by the 
analysis of simpler cases. The mass transport was than 
linked to a two-dimansion.al vertically averaged tidal 
dynamics model for the Providence River. A uniform 
field was simulated, indicating mass conservation 
errors of less than 0.1 %. Modeling coliform 
concentrations in the area revealed a mass conservation 
error of as much as d.5 due to an extrapolated, tima- 
varylnq boundary condition. However, the modal compared 
quite closely to a set of fiald data when no decay was 


specifieci . 

Adiihional effort was devoted to the extension of 
the model to a steady-state application, by replacing 
the time step with an iteration sequence. This was 
verified bv comparison to analytical solutions, and 
demonstrated by application to a river confluence 
situation. Another application of the time-varying 
model was to point sources in Block Island Sound. A 
two-dimensional modal prediction was compared to the 
three-dimensional distribution for the vertically 
well-mixed case, and found very similar after several 
tidal cycles. 
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I. 


INTRODUCTION 


ORIGINAL PAGE IS 
OP POOR QUALITY 


1 


Hater Quality 

"Water quality" generally refers to the ability of 
a body of water to support a healthy natural ecosystem, 
and additionally, to lend itself to the pleasures of 
man. The addition to the water of any substance that 
impairs this ability results in the degradation of water 
quality. Historically, it has been one of the 
"pleasures of man" to use bodies of water to dispose of 
various waste products. Hhen the ocean or coastal 
waters ware used for waste disposal, they were 
considered to be an ideal and limitless sink, capable of 
absorbing the wastes without reducing the water quality. 
Indicators of water quality (dissolved oxygen content, 
biochemical oxygen demand, fecal bacteria, heavy metal 
concentrations) were expected to be reduced to normal 
(ambient) levels due to the great assimilative or 
self-purifying capacity of the sea. 

The growth of population and industrial production, 
especially close to the sea, in the past century has 
begun to strain the assimilative capacity of nearshore 
waters. Many harbors and confined bodies are so fouled 
that the natural aquatic populations are destroyed. 

Some areas that were thought to be in a natural state 
have revealed dangerously high levels of pollutants 


(D*. Some offshore areas close to population centers 
show marked degradation due to offshore dumping (2,3). 
Yet in other areas, similar quantities of waste have 
been dumped with minimal observed ill effects. (4,5). It 
is clear that under some conditions, coastal waters are 
being taxed beyond their capacity; yet if conditions are 
suitable, the water can survive being used as a dump. 
Certainly, the need for v;aste disposal areas is pressing, 
and land disposal presents equally complex problems. 

The critical question, then, is to what extent 
waste disposal in coastal waters can be permitted, 
without rendering the waters unfit for productive and 
pleasurable activity. Three aspects of this question 
must be dealt with: (1). w'hat levels of which 

pollutants will produce what harmful effects? This 
question is the province of scientists, and the answers 
Involve extensive and difficult laboratory and. field 
testing, with results that to date are unsatisfactory. 

(2), ^mat methods can be used to reduce harmful wastes 
to acceptable levels? Here the engineers join the 
scientists to develop or improve methods of waste 
treatment, outfall design, and environmental monitoring 
and modeling, (3), What are "acceptable levels" of 
pollution? This question can only be answered by the 
political process. The public must define, through 

* Numbers in parentheses refer to listed references. 
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government, its economic priorities with regard to 
pollution. The usual practice is to classify waters 
according to their suitability for various activities 
under existing conditions, and to upgrade the waters 
whenever possible (6) . As indicated by the recent 
defeat of sewer bond issues in Rhode Island (7) , the 
public is likely to accept some reduction of water 
guality, rather than bear the cost of cleaning up. 

ORIGINAL PAGE IS 

Water Quality Models OF POOR QUALllIf 

Questions (1) and {?.) , above, form the realm in 
which water quality modeling can be of use. The term 
"model" has various levels of meaning. In a general 
sense, a model is "a conceptual idealization or 
simplified representation of a physical process" (B) . 

This could include physical (hydraulic) models. A 
slight refinement is to define a model as the 
mathematical expression of the physical process; that 
is, the application of the principles of Newtonian fluid 
mechanics with approximations aopropriate to an estuary 
or continental shelf (8) . In the most specific sense, 

"model” refers to a mathematical formulation of the 
physical problem, along with an appropriate solution 
technique. 

Water guality modeling has two purposes; diagnosis 
(identifying and isolating factors affecting water 
guality), and prediction of the effects of changes in 
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the system. A predictive model offers a technique 
widely useful in the design of waste disposal systems. 
The usual objective is to predict the distribution of a 
dissolved constituent or water quality indicator. 
Comparison of the distributions under different 
circumstances can aid in outfall location and design, 
setting water quality standards, choosing optimum levels 
of waste treatment, or anticipating the effects of 
special circumstances such as sewage plant overflow 
caused by storm runoff. 

The mathematical model involves solution of the 
hydrodynamic equations and the parabolic, partial 
differential equation for transport of a dissolved 
constituent, by numerical techniques if necessary. Work 
began on this problem by using a one-dimensional 
approximation, a reasonable simplification for stream 
flow with a plane source of material. This was shown to 
be useful in rivers and in long, narrow estuaries 
(9,10,11). The more difficult two-dimensional 
computations have been applied to estuaries and lakes 
(12,13,14,15,16). In this case, either vertical 
averaging (13,14) or lateral averaging (15,16) was 
employed to eliminate the variation in whichever 
dimension was least important. There remained 
considerable uncertainty about the third dimension, 
which limited the applicability of these models to 
either shallow or narrow estuaries. The solution of the 
three-dimensional system presented a formidable task and 


promised to consume large amounts of computer time and 
storage. 

nany estuaries do not lend themselves well to a 
two-dimensional approximation. Rstuaries with 
insignificant lateral currents and depth variation are 
rare. Although two-dimensional models can incorporate 
those variations, they can shed no light on the effect 
of outfalls at different depths, surface versus 
submerged discharges, or vertical stratification. A 
particular problem posed by stratification is the affect 
of seasonal variation in salinity structure. Also, due 
to the increasing use of relatively unconfined 
continental shelf waters for waste disposal, knowledge 
of the behavior of a dispersant in all three dimensions 
becomes essential. Information concerning vertical 
distributions is of particular importance in continental 
shelf waters: for instance, the settling patterns in 
barge disposal of sewage sludge — a problem of 
increasing concern. 


The Problem and the Method 

The greatest obstacle to developing a three- 
dimensional model is that the vertical velocities are 
largely due to gravitational flow, which can be 
neglected in horizontal flow. Since gravity flow in 
estuaries and seas is mainly due to salinity gradients. 
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it would seem that the hydrodynamic and mass transport 
models must be coupled — a serious complication (8) . 

The hydrodynamic models used to obtain the current 
V( locitv input for the mass transport equation have used 
an approximation of constant density, which is useless 
in obtaining vertical velocities. The approach to be 
taken here, in the absence of a working 

three-dimensional hydrodynamic model, is to neglect all 
yef'^ical velocities, and model the vertical transport of 
constituents by the diffusive mechanism only, which can 
incorporate salinity effects. This will be justified in 
Chapter TII. The model has the capability to utilize 
vertical velocities, however. Models to generate the 
vertical velocity components, currently under 
development, will permit full realization of this 
capability. 

h computer program to solve the three-dimensional 
mass transport aquation by an alternat inq-direction- 
implicit (R.D. T.) finite difference method (17) has been 
written. A tidal current model has been developed for 
the Marraqansett Bay area (IB), based on Leendertse's 
long-wave propagation model (19). This has bean 
employed to supply the required tide heights and 
horizontal current velocities. 

Coliform bacteria concentration, the most widely- 
used indicator of sewage contamination, has been modeled 
only in a very sketchy manner, using a constant decay 
coefficient (20) , or an oversimplified advection model, 

PARE ’• 
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with a constant, two-zone flow approximation (21). In 
addition, a well-distributed coliform data set exists 
for the Providence River estuary at the head of 
Narragansett Bay (22). This area is suitable for 
model d<=velopraent, due to its well-defined physical 
boundaries and sewage sources. This area is therefore 
chosen as the model area, with coliform bacteria as the 
quality indicator. 

At this point, it is well to define the 
relationships of the physical and biochemical parameters 
which form the system to be modeled. The features being 
considered are shown in Figure 1.1, modified from 
Leendertse's model (14). The advection model will be 
discussed in Chapter V; the kinetic model in Chapter VI. 
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Figure 1.1. The environmental model 
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WASS TR&MSPORT EQnATION ANT) SOLUTION 


Basic 'jlquation 

The equation for transport of dissolved or 
suspended matter involves three mechanisms. Material is 
transported by advection due to the mean tidal velocity, 
and by dispersion due to turbulent mixing. Also, matter 
is transported into or out of the field at the 
boundaries, ingested at sources, and removed or 
regenerated by reactions in the field. Using the 
familiar oceanographic coordinate system, in which the 
X- and y-directions are in the mean sea level plane, and 
z is directed downward froti mean sea level, the mass 
transport equation is (23): 



where p is the concentration of dissolved constituent 
u, V, w are the x-, y-, z-directed velocities 
e^, , 0 y, Oj, are the turbulent diffusion coefficients 
in the x, y, and z direction, respectively 
S is the 50urce*sink term. 

The above includes approximations appropriate to 
estuaries and continental shelves. The diffusion terras 


express the time-averaged turbulent flux terras uP^ as 
the product of an eddy diffusion coefficent and a mean 
concentration gradient. They do not include raolecular 
diffusion, which for large bodies of water would be 
several orders of nagnitude smaller than that due to 
turbulence. The velocities u, v, and w are time 
averages over a period much shorter than a tidal cycle, 
but longer than the turbulent fluctuations. 

Spaulding (24) has found it helpful to non- 
dimensionalize the vertical axis. This eliminates 
complicated boundary conditions at the estuary bottom 
where depth variations are considerable, and also 
facilitates inclusion of the tidal variation in the 
overall depth. The parameter used to non-diraensionalize 
is the sum of the mean sea level depth, h(x,y), and the 
tidal height, ^(x,y,t), at any time. The sura is called 
H, A dimensionless vertical scale, rj, is defined: 

•? = _ 2 -? or 2 = h R ♦ f ( 2 . 2 ) 

H 

Thus 17 eguals zero at the free surface and -1 at the 
bottom. A fairly complicated transformation is required 
to express Equation 2.1 in terms of the new vertical 
coordinate. 

Trom Equation 2.2, since H and z are independent, 
dz = H d»^ (2.3) 

Using the symbol S to represent derivatives in the x-y -<7 


free surface 
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coordinates, the relationships between the partial 
derivatives in the x,y,z system and those in the x,y,i? 

system are: 




0 _ S 
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Sx 

H3^ 



(2.J) 

3 - 

S 


S 



(y-fe) 

3y 


^y 

n<y») 





= -i, (^-7) 

H^>7 


To find the relationship betaeen the real Telocity » and 
the dlaenslonless vertical velocity, ID. Eqoatlons 2.4 to 
2.7 are substltated Into the expression for the total 
derivative D / Ot. Since by definition, nf/ Dt = lO. the 
real vertical velocity is related to the dimensionless 
vertical velocity by: 

.au* ^ u + V Ikli) U.9) 

tOHt 

sobstltntlng equations 2.4 through 2.8 into the 
advectlve terns of Equation 2 . 1 . rearranqlnq, and 
multiplying by H, yields; 

jH/o + IpliS. = HT, 

■Jt Sx Jy h 
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where T.^ represents the dispersive terms. 

To express the dispersive terms compactly, let 


Xx 



-Xy Sx H 



Transforming the z-dlsperslve terms by Equation 2.3, and 
the X- and y-terms by Equations 2.10 and 2.11 
respectively, the full equation for mass transport 
becomes ; 



If one Is willing to sacrifice some accuracy, the 
dispersive terms can be expanded and simplified by 
assuming the higher-order terms to be small. The simplified 
form is then: 



Reaction Model Concept 

The source-sink term, S In Equation 2.12 above. 
Includes the reaction kinetics. If any, of the water 




parameter* The composition oT the term will 
depend on whether the constituent Is conservative, 
undergoes simple decay, two-stage, or multi-stage 
reactions. Examples of the source term formulations 
follow: 

(1), Salinity (conservative) (24) 


(2). Conform bacteria (single-stage decay) 


s = Gs + K^c 


(2.H Id) 


(3). D0-30D (two-stage consecutive reaction ) (8) 


DO: S = CgjjQ + Kjj CoBp - 

BOD: S = 0s9oa’ 

where Is the BOD decay rate 
K II is the reaeratlon rate 
Coef Is the DO deficit, Cqq 

Gjoo sources and sinks of oxygen, l.e., at the 
boundaries, 

is waste-load BOD. 


(4), Nitrification (multi-stage consecutive) (8) 
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Organic matter: 

S = CsK, 

+ 



/immonia : ' 

s = 

+ 


a\+ K,^N, 

PI trite : 

s = 

+ 


Kj + Kj 

Pitrate nitrogen: 

s = c,,, 

+ 

K 4 . 

K3,N3 


where is a concentration, and the subscriots 1, 2, 3, 
and 4 refer to or^ranlc, ammonia, nitrite, and nitrate 
nitrogen respectively, 

As the examples above indicate, this problem lends 
Itself to a matrix formulation. /ts suggested by 
Leendertse, this is expressed as: 






(2.15) 


where 3 is the source or sink vector 
is the reaction matrix 
P is the concentration vector. 

For example, in the D0-30D system, this becomes 



[K]P 

4 S 
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Numerical Techniques 


A basic approach to an initial- value problem 
expressed as a set of partial differential equations is 
the finite-difference approximation. This means that 
the differentials are expressed as a small but finite 
step in space, over ■which the valuables of the dependent 
variables are calculates repetitively at similarly small 
time steps. The followinp: integer subscripts will be used 
to denote spatial Increments and the time step in the 
corresponding directions; 

rr. ; >:-directlon 
k; y-directlon 
n: f^-dlrection 
1 ; time 

?or parabolic, and also for elliptic, equations, 
the alternating direction method (1?) has been 
effective. This involves splitting the equation into a 
different level for each direction, and advancing in 
fractions of time steps. The advantages are 
unconditional stability, suitability for an implicit 
solution, and the ability to solve by a tridlagonal 
matrix technique. This method was extended by Douglas 
(17) and Douglas and Gunn (25) for three-dimensional 


computatians, ani will be employed in the present study. 
An additional refinement is to choose a space-staggered 
grid system, as shown in Figure 2.1 (17). The purpose 

of the staggered grid is to have each variable centered 
between the ones upon which its calculation most 
depends. The distances are expressed as the products of 
the aforementioned integers, m. It, n, and the spatial 
grid sizes Ax, ^y, The mass density, p, is now 

called P, to indicate that it is an average value over 
the grid volume, rather than a point function. The 
functional relation is now exnressed as (24): 

P(x,V,«7,t) = P(mAx,ltAy, n47,14t) = (7*1*5) 

The following difference operators are defined: 

II L 

Hn,K,n = ’PyY,,k,vn+^ 

The first equation below calculates the x-directed 
concentration variations from time-level ji to 
These values are then used to calculate the y-directed 
variations from time-level Ji to Finally, the 

) 7 -variations are calculated from level 1*3 to 1*^, 
using the previously obtained values. This completes 
the time step, rising the approach of bouglas (1 ) , the 
finite-difference form of Fquation 2.12 is the following 
three-equation system (24) : 


h“ Mean sea level depth 
tidal height 

u, ▼, Hi- velocities in x, y, directions respectively 


P- 

^X' 


mass density 

D , D^- Dispersion coefficients 


Figure 2.2. 


Three-dimensional space-staggered grid system 
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2 2 - 
A “t ^ 

+i5,(DxH<S/P)/-<^y(VHp/ - (P>H'’')^ +<5v,(DyH^v('’)) 

.<S',C^^.?CP))^+SH^ (2.17) 


(^p/ r= _i ^y(UHP)^ ^ 

.1 4(Ih<S,(P)/- i ('/"P)-"'’ -i^y(''«P)^ ^^(OyH^vfP^^'’ 

5y(iV H (u)HP)^ + <r,, (^ (p))^ (1.18) 


(HP^- (WP)^ = (uHp/\ JL (5,(uHP)^4iX/D,H5y(P))'^'^ 


2 

j. 
2 


c^xCD^H5,CP)y -4 CVHp/-^-l ^y(VHP)^+ii^VW) 




i4 

£^-1 


A simplification can be made bv subtracting 
KquatiM 2.17 from 2.18, 7n« 2.18 trq. 2.18. This 
yields (24): 

X((ytp)^’i(np/j =-^xCuHP)'^^-5«CuHP)-‘ ^ 

rlcDyH^,(P)/^7 (D,H5.(P)/ -2 5y(VHP)^-2 J.(p> ^P) 
+l^y(t>yH^y(P)/-l2<5''?(^ ^ 

2. ((HP)^*^4Mp/'^) = -5y(VHP)^"^-^y(VHP)^ 

+ .Jy (Dy H Jy (P))2-*§ _ <Jy (DyHi'y(P))^ 


(2.21) 
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2- + t5''7 (ujHP)-^ 

Th 9 source-sink ter«, (HS) In Siimtion 2.12, is 
appr-.«i»atea by the cenction latrtx schene. In finite 
difference notation, this tern becones 

'vi+4 


M OjiAK ii I (2.23) 

^ / . 1 1/ r%X X-+1 


(HS;)'^'^= r (HK-.j + (HSiV 

j = ' 


■ her.-- the subscripts i and 1 distinguish elenents of the 

reaction matrix. ^ 

Tn the above, the value of ?f is unknonn. If the 

previous value of P is used, accuracy is 
alternative is to estimate a new value for p/* from 
Equations 2.20, 2.21, and 2.22. Then the solution must 
be readvanced in time using the e.stimate. This restores 
the accuracy, but doubles the connu tationa 1 time. The 
additional accuracy gained bv this estimation normally 
does not iustify the additional computational time 
requirement, and is therefore not u.sed in this study. 
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Qualitative definitions of diffusion and dispersion 
have been suggested by Holley z "let diffusion 

refer to transport in a given direction due to the 
difference between the true convection in that direction 
and the time average of the convection in that direction 
Let dispersion refer to th- transport in a given 
direction due to the difference between the true 
convecfion In thnt Jlcontlon anrt t.h- snatUl a.etaqn of 
thn oonveotlon In that direction.” FirinflV, Itffuston 
is <iiin to nolnonlar and torbnlsnl notion, and dispacsion 
is duo to tho variation of the nsan velocity across a 


section . 

A diffusion coefficient is given by 


n = J. d ( ) 


(3. 1) 


2 dt 


where cr^ is the mean square dispersion of particles 
(27): in other words, the standard deviation of the 
concentration distribution in the dye plume. This 
consists of a molecular and an eddy diffusion term. The 
former is negligible, esDscially in ooen waters, where 
the molecular effect is several orders of magnitude 
smaller than the eddy effect. 
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Prandtl expressed the diffusion coefficient in 
t.erins of fluid turbulence (27): 

where 1>< is a mixing length, and u’* is the mean square 
velocitv fluctuation. The mixing length can be defined 
in sev-ral ways, and will be discussed in a later 
section. The velocity fluctuation is the difference 
between the velocity at any instant and the mean 
velocit V". 

( 3 . 3 ) 

u = u + u’ 

The mean square dispersion is aiven at any time by 
the concentration of the dispersant at a distance x 

• , f"" -z j 

^ ^ ‘•x ri.41 

where b^ is the initial width of the natch of 
dispersant, and c^ is the initial average concentration 
therein. Foxworthy found that as the time scale of the 
process increases, the value of cr is time-dependent 
(29). For small, intermediate, and large time scales, 
was found proportional to t^ , t^ r and t respectively. 
Thus the diffusivity is determined by different 
functions of time, until the time scale is such that the 
patch has become larger than the characteristic eddies. 
This is called the asymptotic nhase. 


Oceanic Diffusion”" Horizontal 
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The horizontal eddv diffusion coefficient. 


in the 


taken as a 
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absence of a stream flow, is generally taken as a 
function only of a length scale, which increases as the 
patch disperses. The simplest definition of the scale, 
1^, is the width of the patch at a given time (27). as 
a mixing length, is defined as the characteristic 
distance traveled by the eddies before losing their 
identity. This, however, is not a very useful 
definition. Talcing the scale as the patch width, this 
is variously defined according to the dispersion . 

For lateral diffusion (27) : 


ly= 4cr 


(i.'j) 


or, tor a line source, so that ly = b at x =0: 


2VT cr {^*6) 

For radial diffusion, Okubo gives (30): 

3<r., (3.7) 

where the subscript rc signifies the radial 
distribution. The above all specify that 05* of the 
diffusant is contained within a distance of 4g , , 

or icr^j.. 

Oceanic turbulent diffusion is considered to 
consist of three phases, corresponding to the 
aforementioned time scales (29) : 

(a) , eddies larger than the initial dye patch, 

(b) . eddies comparable in size to the dye patch, 

(c) . eddies smaller than the dye patch. 

(c) is the asymptotic phase, and produces the greatest 
rate of mixing. In finite difference modeling, fhe 
requirement of an initial average concentration 
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throughout a large finite volume makes it reasonable to 
assume that the diffusion is in the asymptotic phase. 

If the model scale is such that eddies larger than a 
single grid exist, the eddies then appeals part of the 
velocity fieli. As the resolution of the hydrodynamic 
model is improved, the dispersion terms become less 
important. 

The basic relation considered to hold for oceanic 
turbulent diffusion, in the asvmptotic phase, comes from 
Richardson's ]aw for atmospheric diffusion (31): 

where is the horizontal eddy diffusion coefficient, 
and is the dissipation parameter. This is called the 
"four-thirds law", Many field measurements of the 
diffusion coefficient have been made. The general 
procedure is to use a fluorometcr to measure the 
concentration of dye tracer at a known distance and time 
from its injection. The coefficient is calculated from 

- ^i'^) (3.9) 

(■^2 ” ' t ,) 

In the attempt to determine the validity of the 
four-thirds law for diffusion in ocean waters, and to 
determine the values of the coefficient, compilations of 
the field data have been made by Yudelson (32) and Okubo 
(30) . Their diagrams look very much the same. But 
Oku bo, in the mors recent report, has the more 


(26) : 


Dv = ^ 


I 
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enliqhtoning c:onc:lusions. 

Figure .1.1 shows rlif f usivitv versus length scale 
for data obtained in 44 surveys (12) • The locations 
include lakes, estuaries, bays, coastal waters, shallow 
seas, and deep ocean. The lines drawn are the 
four-thirds law for different values ot the coefficient 
(0.1 to O.oni cmVsoc.). The scatter spans two orders 
of magnitude. This is not entirely surprising, 
considering the variety of locations. However, values 
for different ocean locations are as widely scattered as 
any. yudelson notes that many points fall on the 4 /I 
power lines for otjj's of 0.02, 0.01, and O.OOS cm. /sec. 
This happens in th° middle range of 1^, from 10m. to 
10km., which is the most useful range for modeling. 

Okubo (10) shows a diagram. Figure 1.2, very 
similar in appearance to that of Yudelson, The data are 
from twenty investigations, all slightly more recent 
than those cited by yudelson. between the two lines 
delineating the four-thirds law, a shift from left to 
right with increasing scale is apparent, indicating that 
the overall exponent is le.ss than four-thirds. This 
trend is also apparent, although not noted, in 
Yiidelson's diagram. 

Fitting ^ine to his points by eye, Okubo offers 

(30); 

P^= 0.0101 Ix’*^ (d*10) 

which has also been drawn in on Yudelson 's diagram. 

Okubo suggests that the four-thirds power law is valid 
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only locally for some lenfcth scales. This is shown by 
dotted lines fitted by eye, v;here: 

0.008 for the range 1 . 5 x 10 ^ cm 1 / 2 x 10 ^ cm, 

0^0 0.01 for the range lO^^cm, < ly< 1 . 5 xl 0 ^cm. 

0,03 for the range lO^cm, < ly<l 0 ^cm. 

Obviously, the data sets used by Yudelson and Okubo do 
not coincide very well (Okubo*s line, Equation 3.10, is 
to the right of almost all of Yudelson*s points). In 
the Interest of precise comparison, both data sets 
(ref, 32 , pp.k -2 to iU 8 , and ref 30 , pp ,793 to 795 ) 
have been subjected to least-squares fitting, and are 
dravm in ?lgure 3 . 3 , Also shown is a four-thirds law used 
by Koh and Chang (33) in a barge-dumping model, in which 
0 ( 5 = 0.00015* ?or the data of Yudelson and Okubo, only 
the points in the useful range for modeling, 1 meter to 
100 kilometers, were used. The equations obtained, and 
the standard deviations in orders of magnitude (since it 
was a lo.'-arithmic operation), are as follows; 

Okubo : 

1.0384 

0^ = 0,0366 1^ , lyin cm., standard deviation 0,051 

0.00136 ly , 1^ in feet (3.11) 

Yudelson: 

I /444 

iJj, =0,0655 l^’ » ly 

0.00333 ly^^ , lx in feet (3*12) 

standard deviation 0,3056, 1^3 data points 
Composite; 

Oy= 0.18 1^ , ly in cm. 





D^= 0.00S9 1^ , f«et 
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(3. 1.3) 

st^ndarc! deviation .3445, 198 data points 

Equation 3.11 is sliqhtly different from Equation 
3.10, due to the dropping of data for greater than 10 
cm. The slope of the composite line is reduced because 
Yudelson's data included more points in the lower range 
of I than did that of Okubo. Perhaps because of this, 
Okubo's relation is preferable for offshore models. 

Okubo was very demanding in his selection of data, 
requiring dye release to resemble a point source, and 
the distribution to be measured according to certain 
specifications. It is not likely, though, that this 
justifies ignoring Yudelson*s data. The composite 
relation is probably most representative, although 
Okubo’s is more attractive. If there were agreement 
about the localities in which the four-thirds power law 
holds for a certain value of this law might be 
useful. But that does not appear to be the case. 
Equation 3.13 indicates that may actually be directly 
proportional to 1^^. 


Studies of dispersion in unconstrained waters 
indicate the variation of the dispersion coaffioient 
only as a function of the scale. In rivers and 
estuaries, other factors are included: mean flow 

velocity, the bottom roughness as expressed by the Chazy 
coefficient, the depth of the water, wave action, 
hydraulic radius, and tidal period. An approach to 
modeling horizontal diffusion in open waters would be to 
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assign a value of to ^ach spatial grid by one of the 
equations 3.11, 3.12, or 3.13. The lenqth scale would 
be twice the distance from the grid to the point of 
origin of the pollutant. This means that dispersants 
from several sources would have to be modeled as 
separate constituents, using a different set of 
dispersion coefficients for each source. 

vertical Diffusion coefficients 

Diffusion in the vertical direction is restrained 
by stability, by the bottom and the free surface, and by 
the magnitude of the vertical component of flow, which 
is normally very small relative to horizontal flow. 

Thus it has an effect several orders of magnitude 
smaller than that of horizontal dispersion. The latter, 
however, is still a small effect when dispersive 
transport is compared with advective transport (14). 

This is not likely to he the case in the vertical 
direction, since the vertical flow is also very small, 
vertical diffusion and vertical advection are verv 
closely related affects, since they are driven by the 
same force (instability). Since diffusion is also 
highly dependent on wave action, it is likely to be more 
important than advection. This indicates that ignorance 
of the actual vertical advection- or the inability to 
simplify it to finite-difference grids, may be partially 
overcome by using diffusion alone to model vertical 
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transport. 

Th?? simplest method of calculating the vertical 
diffusion coefficient, D , is correlation with the 
Richardson number t 

This incorporates stratification and the rate of lateral 
flow. However, since vertical diffusion is a sensitive 
procrjss and difficult to measure, it has proved very 
difficult to define this relationship. All the 
correlations listed in Table 3.1 have been offered 
(33,34). Calculation from the Richardson number would 
require salinity and temperature, or density profiles 
from the region of interest. 

The only general aoproach to determining T>^ iii to 
measure the vertical distribution of some constituent in 
thve region of interest, and note the corresponding 
physical factors, such as surface waves, the salinity 
profile, depth, and currents. This has been done for 
heat and for suspended matter by Sastry and Okubo (35), 
and Tchiye et al. (36), respectively. 

Ichiye et al. obtained an in-depth picture by 
measuring the distribution of suspended matter in deep 
waters, for nine widely separated stations in the 
Carribean. Yet not enough data were obtained to relate 
eddy diffusion to static stability conclusively, was 

found to vary both horizontally and vertically. In the 
upper layer of the sea, to 150 m., 0^ varied from 0.7 to 


Table 3,1. Correlations of J x^ith Richardson I'jucber 

or Density Gradient 6 

^zo = 3z R; = 0, the neutral case 


p is f- proportionality 

constant 


P.ossby & Montgomery, 1935 (33) 

(3.15a) 

Rossby & Montgomery, 1935 (33) Dz= Oz®(l 

(3.15b) 

Holzmt.n, 19^3 (33) 


(3.15c) 

Yamamoto, 1959 (33) 


(3.15d ) 

Mamayev, 1958 (33) 


<oJ 

• 

CD 

Kunk 1; /xnderson, 1948 

(33) 

/3=3.33 

(3.l5f) 

Harreroes, 1968 (33) 

-3 -"i -j 

Dz •=* S'y 10 £ cwi /Stfc. 

5ylO’V € <• ISy I0"^cv>t"‘ 

(3.15=^) 

Kolesnikov, 1961 (33) 

Dz = ^^Z6 % J 

/3=8.3yl0'^ 

Dzo = 2 (3 = lO.Oy I0-5' 

(3.15h) 

Koh and 3'an, 1969 (33) 

4ylO"^i 6 ^ \ 0 '^ 

(3.151) 

Guttman and Huang (34) 

D_ =uLR;"^^2, 

(3.15J) 


L ^ lciA.<^+U scale 

Uj vcloc'i4v| scale a"f wo't'iow 


8 cm^/sec. 

13 cm^/sec., 
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In intermediate depths it was between 5 "^nd 
and from 0.2 to 2 cmVsec. within 80 m. of 


the bottom. 

liullenberpr (37) offers a definition of and a 
simple relation based on measurements in shallow v;ater, 
l4 to 30 meters. He finds that stratii ication and •rind 
velocities have marked influence on vertical diffusion. 
The survey was based on the following; definition: 


= 


—dill— l«f4L) 


(3.K.) 


where h^ is the half- thickness of the dye layer, above 
or below the point of injection. Kullenberg derives the 
following relation from his experiments: 


D, ^ 8.^ ylO'® 

^ “ dz. 


( 3 . 17) 


where v; is wind velocity, m /sec. 

q is the horizontal velocity vector, m /sec. 

.. is the stratification number, t 1 - 

3 pAz sec®- 

3ome sort of vertical velocity profile must be 

obtained to estimate dq /dz. The obvious criticism 

Of this relation is that it implies no vertical 

diffusion in the absence of wind, and infinite diffusion 

in unstratified waters. The variation of vrith depth 

is deuendent on the functions chosen for 11^ and dq ,/dz. 

bastry and Olcubo (35) suggest a method for 

obtaining heat diffusion coefficients, which could be 
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applie.^ to concentration by analoqy. The metho-l 
requires knowleiqs of the temperature (or concentration) 
profile, anfl the absolute vertical velocity profile, in 
the waters of interest. The qoverning assumption is 
that the sura of vertical advection and vertical eddy 
diffusion is the same at any depth. A value of D, is 
estimated for a certain depth, say 100 ra. Then: 

D,-_(|Ey' [.T-wT,„.(D.g)J 


Tf the temperature T and the vertical velocity w 
are known at all depths, a profile of can be 
obtained. The assumed value of (D^),oo 

on whatever information is available, but values which 
would make 0^ negative at some depth, or unreasonably 
high near the surface, compared to measurements found in 

the literature (33). can he eliminated. 

The most flexible and most comprehensive method is 
that offered by Pritchard (38). Based on measurements 
in the dames Piver. it calculates 0 from the Richardson 
number, wave parameters, horirontal flow, depth, and 
empirical coefficients. Although the coefficients have 
only been obtained for the James River, the formulation 
was designed for partially-mtx^d estuaries, and should 
be fairly general. Since field measurements of 
stratification are required for all methods, some 
diffusion measurements taken at the same time would 
permit adjustment of these coefficients. Pritchard's 

formulation is: 


Pz = IpmUh-z)^ R;)'^ (3.19)^^ 

where v]f , , «ip ara aiiustable coefficients 

7 is distance downward from the surface 
WH is wave height 
NL is wave length 
WT is wave period. 

This formulation has been used by Spaulding (11) 
without changing Pritchard*s coefficients. No matter 
which method is used to calculate , some field 
information is necessary. In the open ocean case, if 
diffusion measurements cannot be made, values must be 
assumed based on the most similar cases in which *’ 
measurements have been made (31) . 

Horizontal Dispersion in Estuaries 

Mass transport in estuaries is usually dominated by 
currents caused by tides or river inflow, and the effect 
of bottom and shore is more pronounced than in offshore 
areas. For these reasons, a length scale is not 
considered sufficient for calculation of and Dy . 
Methods that have been used for estuaries are based on 
the mean current velocity, and are therefore said to 
calculate dispersion coefficients rather than diffusion 
(14,39) . 

Holley, Harleman, and Fischer (39) present a rather 
complicated scheme for tidal estuaries. It appears to 
be applicable to wide coastal bays if the hydraulic 


radius is takan as equal to the depth. New physical 
factors taken into account here are the oscillation of 
the tidal flow, the Chezy coefficient, and the width of 
the body of water. 

The method begins with Elder’s empirical relations 
for longitudinal, lateral, and vertical dispersion (40); 


Px = 5.93 h,u* 
Py= 0.?3 hjU* 
nj,= 0.067 lyi* 


(3.20) 

(3.21) 
(3. 22) 


Where 




^ 1 

•ni 

u* = 

^ u 

■ K 

u I 


^v. 


is the Chezy coefficient 
hj,is the hydraulic radius, which is taken 
equal to the depth. 

Two dimensionless variables incorporating the tidal 
period are introduced: 

= X./J^ (3.23) 

/ Dy 




(3.24) 


where h is the estuary half-width 
h is the depth 
T^ is the tidal period. 

Holley et al, then produce a formula relating the 
diffusion coefficients in each direction to the 
coefficent for an infinitely long tidal cycle, n< 50 . 
Fortunately, for vertical variations in the mean flow. 




thin reduces to: 


5.93 ha* 


( J. 25) 


I 

m 

1 


Dy - Dooy ~ 

For transverse variations, 

■= 10 (T„')^ 

Doot (3.26) 


1 

I 

I 

I 

i 

1 

1 

I 

I 

I 


D-t = 

(-3.27) 

33 where u"^ is the cross-sectional average of the 
sq aared velocity f luc t uat ions, and 

u'* = 0 ( 5 ( 2 - :^Z) sin (27[t) (3.28) 

Very little is actually known about the likely values of 
u'*, and a value of the dimensional constant must be 
assumed. In the above, and represent the 
longitudinal dispersion due to vertical and transverse 
velocity variations, respectively. If the velocity 
distribution is known, Dy can be taken equal to the 
larger of or . Tf not, it is recommended that 
be taken equal to D^. 

Leendertse (14) offers the most appealing method 
for calculating dispersion in an estuary. (!e also 
starts from Flder's equations, 3.20 and 3.21. The 
trouble with these is that they indicate no dispersion 
in the absence of a mean flow. However, it is clear 
that dispersion in a natural body of water will still be 
considerable, due to turbulent motion and wind effects. 
Leendertse assumed the following functional 
relationships: 

0x= ♦ ’’w 
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+ n 


w 


(3.30) 


«horP n ir, a diffusion coefficient dependent on wa ve 
and .ini conditions, and on the lateral diffusion. 

For use in this development, which is to be applied 
to a partiaily-mixed estuary, Pritchard’s method for the 
vertical diffusion coefficient will be used (Fquition 
3.19). ^he horizontal disp sion coefficients can best 
be calculated from Fouations 3.29 and 3.30, using 
Equations 3.20 and 3.21 as the functions f, and f ^ . fhe 
value can he taken as a diffusion coefficient, as 
calculated from Equation 3.13, using a typical width of 
the estuarv as a length scale. This would oermit 
adaptation of the same method to offshore areas, whether 
or not the mean flow is dominant over turbulence. 

Table 3.2 summarizes the aspects of the different 
methods presented. The information listed must either 
be gathered in the field, or assumed similar to data 
alreadv available (if it is not self-evident, as would 
be values for the length scale) . 


ORIGINAL PAGE IS 
OF POOR QUALFIY 


TABLF 3.2. Computational methods for diffusion 
coefficients 

Reference Equations Required Information 


Horizontal : 

1. Power laws (27,26,27) 

11, 12, 13 

2. Elder (39) 20, 21 

3. Holley et al. (39) 20-28 

4. I.eendertse (14) 29, 30 

Vertical : 

1. Kullenberg (37) 17 

2. Sastry and Otubo (351 

18 


3. Holley et al. (39) 25 

4. nichardson no. (33,30) 


. "ritchard (38) 


15 


Ik 

h , U , 

h, u, Cyj, b, T, u* 

h, u, 

M, dq /dz 

vertical velocity 
concent rations 
at some depth 

h, u, 

E‘, f Ojo , 6, ^ 

R. , u , Z, WH, WT , WL 
^ p / ^ 9 ) ^p 


5 


19 


rv 


TSVESTIGRTTON 0? COMPUTATIONAL ASPECTS 


Th?> f inihff-di fferenca equations presented in 


Chapter II are, of course, only approximations to the 


mass transport equation. It is essential to Icnow the 


affects of the approximations. In or<?er for the model 


to be useful, the numerical and true solutions must 


converge as the grid size and time step are decreased 


For simplified cases, analytical proofs of convergence 


and stability can be made. For the more complicated 


system under consideration, it is desirable to 


demonstrate the computational veracity with some 


simulations. This has not been done previously because 


of limitations in computer size and speed. However, the 


time vas taken here to attempt verification of some 


predicted computational effects. 


Convergence 


In order to demonstrate convergence of the 


numerical solution with an analytical solution, a 


rectangular basin with a constant, unidirectional 


current and an instantaneous plane source of material 


was simulated. The solution for these conditions is 


given by Diachishin (41), as follows; 




a 






l-N 


Ml' 
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P= ^'.1 evp L 4Pvi -1 ' 

-1 7T Dj( t 

where gp| is the instantaneous plane source strength 
per unit area, 

lijj is the longitudinal dispersion coefficient, 
y is the downstream distance, 
ri is the constant stream velocity, and 
t is the time elapsed after injection. 

The plane source was simulated by a row one grid wide, 
across the entire width and depth of the basin. Into 
this row a relacively dilute constituent was injected 
over one time step, so that the total mass injected per 
cross-sectional square foot was equal to . The 
distribution was printed at chosen time steps as the 
constituent was advected and dispersed downstream. It 
is expected that as the time step and grid spacing are 
decreased, the finite-difference and analytical 
solutions will converge. Thus the simulation was 
performed for both large and small time steps and grid 

spacings. 

The results of the simulation are shown in Figures 
4.1 and 4.2. Tt is obvious that the simulation is very 
poor for the large values of At and dx, and that the 
solutions are converging to an acceptable representation 
as the sizes are reduced. One may note that even in the 
better case, the simulation is still rather poor after 
the shorter simulation time. This is because the 
material is already distributed over one grid Ax in 
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t=750 sec. 


t=3250 sec 


Analytical solution 
. Sinmlation 
gj^= 12500/ft2 
Ax = 1000 ft. 

At = 250 sec. 


Figure 4.1. Plane-source Simulation. P(x,t) vs. x for large time, space increments 






at the source 


hs the 


width, rather than being planar, 
constituent proceeds downstream, this error becomes less 
signiticant. 


Stability and Accuracy 


Stability is the requirement that errors introduced 
in the comoutational method do not amnlify in an 
unlimited manner (19) . A definition of stability is as 


follows (15): 

« 

T f n ^ 

di f ference 


solu t ion , 


is the theoretical solution of the finite 
equation, and ^ is the numerical 

•A j V\ ^ V| 

than stability exists if the difference 


(error) , 


~7 


(4.?) 


remains bounded as 2 increases. 

For a parabolic or elliptic linear equation with 
constant coefficients, the alternating-d irection- 
implicit method has been shown to be unconditionally 
stable (17,42,41). Spaulding (15) has proved 
unconditional stability of his method for the 
two-dimensional laterally-averaged mass transport 
equation, with the restriction that the dispersion 
coefficients are constant. 

According to the von Neumann method of analysis, 
the error (Equation 4.2) is hra rmonically 

decomposed into the error function 
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E W - i Aj 


(4.3) 


wherp the erequenc:y /3j is arbitrary, 

A] is the amplitude, 
i - square root of -1. 

only a single co.ponent, 1 = s, neea be dealt with. to 
follew an error an tine increases, the error being zero 
at tine t = 0. the solntion of the finite-difference 
Gouation is talcan as 


The e 


irror not grow with tim- if the criterion 


( I - * 


(4.4) 


is met. For the restricted cases described, expressions 
for the para.eter e*'^ can be obtained, and can be 
oonpnted to meet this criterion, thereby proying 
uncon ditional stability . 

Although the equations of Chapter II do not meet 
the above restrictions, and stabilitv has not been 
proved for this case, there is numerical evidence that 
stabilitv arguments will remain valid (17,41). In 
applications to be described in Chapter VIT, the 
computational method exhibited stable behavior at all 

times, when correctly posed. 

in a simpler snnsn, r-tability In nlso go.nrnad by 
thn dinopsioplenn paraaetnr n^ , which is considered a 
basic stabllUy indicator of finite dlfferenca methods 


The criterion is: 

\ ^\ - ^ (». 5 ) 

where p* .lepends an the type of differencing system 
employed, and is typically of the order of unity. This 
normally requires that a particle must not be advectel 
across an entire grid during ono time step,^nd thereby 
lost. A similar, but lass demanding criterion involves 
the dispersion coafficiant: 

^ F (4.6) 


Dx 

(A 

which is less demanding because it is easier to meet for 


real bodies of water. 

The accuracy of the alternating-direction implicit 
method has been calculated by Douglas (17) to be of the 
order of (Ax^+ ) , where Ax is expressed as the 

fraction of the model length spanned bv one grid, and At 
as th<=‘ fraction of the simulated time spanned by one 
step. The desired accuracy, combined with tha criteria 
of Equations 4.*i and 4.6, will give an idea of the time 
and grid increments required. However, the geometry of 
the water body may impose stricter limits, and, as 
described below, there are other numerical effects which 
require consideration. 


Dispersive and Dissipative Effects 


Stability and convergence alone are not sufficient 


conditions for a useable modal. The additional 
numerical effects of concern are dispersive and 





dissioative chara:: teristics. These are defined in terms 
of a number of superimposed Fourier series, of different 
frequencies, which constitute the spatial variation of 
mass densitv. The dispersive effect is displayed whan 
the components in the computational model propagate at a 
speed different from that of the analytical solution. 

The dissipative effect is evident when the components 
decay in amplitude without any physical reason. Since 
decay for physical reasons, such as dispersion and 
biochemical reactions, will exist in the real problem, 
it is essential that the dissipative error not be large 
enough to be confused with any real decay. 

For this purpose, analytical investigations of 
these effects have been performed (14, IS). It was 
necessary to analyze a simple one-dimensional mass 
transport equation for a rectangular basin. Therefore, 
it is desirable to compare the performance of the 
proposed finite- difference equations to that predicted 
for the simpler case. This is done by simulating a 
concentration wave in a rectangular basin for a number 
of different values of the dimensionless parameters 
and Pv ^t - sufficient to construct figures similar to 
those given by Spaulding (15) and Leendertse (14) . 

The mass transport equation that was analyzed is: 


+ « ^ - D, 


= 0 


(4.7) 


where u and are constant 


The analysis of Leendertse 
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(14), which is samraari29(^ here, uses ^ nmlhi-oparation 

ji-t-l £.+2 

method, obtaining equations for and . The 

solution in terras of a Fourier series is: 


PCv.-t) = I Pf expUi'TiX +u)jl)] 

j 


(4.B) 


where (/)i is the frequency of the ith component, 

(31 is the wave number, 2n/ T., 
is the complex amplitude. 

The finite difference equations are abbreviated as 


-* 1+1 
p 

Xi P**" 

'w\ — 

# t+1 
>^2 "Pwi 


(4.9) 

where 

X, = ^ 

I t L A + B 


(4. 10) 

and 

A = 51 in (cr A x) 

AX 

B = ^ 5iV\^ 

(^yf 

(4.11) 

obtained from the full set of difference equations. 

The quantities needed to indicate dispersive and 
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dissipative effepts are the .olalus ef the prepa,atlor, 
faster and the phase shift. The propagation faster is: 


i-r ( ka) = 

' ^ AY ' 


exp L i t J 
exp [ I (u)-t + x) ] 


(4.12) 


ahere «• is the frequency of the coaputed «a»e, and up is 
the frequency of the prototype .aye, after the wave has 
propauated one wavelength. The rate of propagation of 
the real wave is <ruat: the coaputed speed Is the real 
part of Equation a. 10. Thus for the two operations of 

. • 4 -jao of compufed to theorotic-^l- 

the c'>mPUtation, tha catio-v t: 

wave speed are: 


R = •Vav\ V T 

CT U 




I 




(T U 


ind the total ratio i; 

R - i Rz) 


Th“ a 


nalytical r,olution for the araolituda of the 


,ass density fluctuation is 
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p(y,0 = exp (4.15) 

Th= modulus of the propaqation factor is found to 

be: 



exp Dx 


H 

z 


(4. 16) 

where n is the number of operations used while the wave 
prooaqates one wavelength. 

To determine the dispersive and dissipative effects 
of the method of solution given in Chapter II, the 
initial mass field was set to a sine wave variation in 
the x-direction, which was renewed at the source as it 
propagated downstream. After it had propagated one 
wavelength, the amplitude and phase shift were noted. 

The modulus was calculated from Equation 4,15 and the 
observed amplitude. By close comparison of Figures 4.3, 
4.4, and 4.5 with those prepared from the analysis (15, 
14), it is evident that the full set of equations 
performs very much as predicted by the analysis of the 
simnier equations, 

Leendertse's observation that there should be at 
least ten grid points per wavelength, to insure that 
dispersive and dissipative effects are negligible, is 
borne out by this study. When the wavelength of concern 
is the tidal wavelength, which is of the order of 200 
nautical miles in Narragansett Bay ( 18 ) , this condition 






Figure 4.4. 


Dissipative Effect. Modulus vs, 

ZT 

Modulus, |T ( W)| = Amplitude, model 

I ^ I Amplitude, analytical 


after one period ^ 

u 
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is easily met, for a total model length shorter than the 
wavelength. A couple of further observations are in 
order, with regard to the behavior of the computed wave. 
First, what has been called the dissipative effect turns 
out to be a divergence, that is, the failure of the 
amplitude to decay as fast as it should, according to 
Equation 4.15. This is in accordance with the behavior 
predicted by Equation 4.16. second, the phase shift was 
observed to occur as the wave was generated at the head 
of the model basin. Only the wave at the head was 
shortened; waves downstream lagged in time but did not 
shorten or propagate more slowly than the current field. 
This is expected, as the number of waves must be 
conserved. These properties are illustrated in Figure 

4.6. 

ORIGINAL RAGE IS 

Discontinuities pQOj^ QUALifY 

A’hen a discharge occurs in a stream flow, a 
discontinuity of mass density occurs uustream from the 
source. Some matter should be taken upstream by 
dispersion, but the computational method is unable to 
properly represent the discontinuity. This is because 
its Fourier decomposition consists of waves which are 
too short to be encompassed in the grids, and thus a 
disturbance is generated at the source. The effect of 
this is to underestimate the influence of dispersion, 
and negative values of mass density may be produced, as 





i 



«: M 


can be seen in Figure 4.1. 

When the location of the discontinuity is known, 
upstream flux differencing, as described by Leendertse 
(14), can be used. This will increase the dispersion 
enough to siipress the disturbance. However, for a 
three-dimensional problem, a more general approach is 
needed. There is the problem of a polluted stream 
entering a larger hay, which can produce a line of 
discontinuities. Also, reversal of the current field 
after slack tide can produce the same kind of 
disturbance. To handle these conditions, Leendertse's 
method of adding artificial dispersion at extreme 
concentration gradients can be applied in all three 
dimensions. The previously calculated dispersion 
coefficients are adiusted as follows: 




where is an empirical coefficient. Similar 
expressions are used in the other directions. It can be 
seen that the amount of dispersion added depends on the 
mass density difference between the two grids. The use 
of this method adds to the computational time, but 
reduces the generation of negative densities. If the 
scheme is not used carefully, though, it has the effect 
of flattening real peaks by adding too much dispersion. 
Therefore a numerical study was made to determine an 




Again a plane 
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'•optirmim" value of. the coefficient e^. Again a plane 
source in a rectangular basin was simulated. The 
coefficient e^ was varied while other parameters were 
held constant. The results are shown in Figure 4.7. 

The optimum value of ej appears to be about 0.20. 
The effect of flattening the peak at large values of ei , 
and the upstream negatives at smaller values, are 
clearly visible. Tn the real-world simulations 
described in Chapter VTT, it was found that values 
larger than 0.20 could be used without flattening the 
peaks excessively. This is probably because, with 
dispersion acting in all three directions, the 
concentration gradients and hence the amount of 
adjustment made are reduced relative to the rectangular 

basin case. 


J 





V. 


TWa-ni!iBNSrOMAL VERTICAI.LY-nVEPArt.ET) TIDAL MODEL 


Th'^ most usaful mathoil to date tor the computation 
of hydrodynamic input for the water quality modal is a 
two-dimensional, vert icallv-averaqed tidal hydraulics 
model. The method was developed by Leendertse (19), and 
applied by Hess and White (18) to Narraqansett Bay. Tha 
program was applied to the Providence River area to 
determine the circulation and the tide height 
information. The use of this kind of model raguires 
that the actual vertical variations of the velocities be 
small (19). 


Equations and Solution Method 


The system of equations to be solved consists of 
the Bulerian Navier-S takes momentum equations, and mass 
conservation for incompressible flow, as follows (18): 


9iA , 
d-b d X 


a ■- - I + -f V + — f a 

d\f 9z pdx P dy Dz. ' 

(ff.i) 



Duv 

17 


4 Jwv _ _ i 
d y cTz P Oy 
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d\Aj , ^ Dww/ - — I — q + _L f cLiiS + + 3 'Tt zN 

Ot c)x c?y p \ c)x d\j Dt. J 

(5.3) 

^ -V ^ = 0 (5. 4) 

3y d-z. 


Tht? horizontal velocities are averaged by 
integrating in the z-direction, from -h at the bottom to ^ 
at the surface. The z-momentum equation is reduced to 
the hydrostatic equation, making the 

Boussinesq assumption that pressure varies only with 
depth. The pressure at the surface is assumed constant. 

Bottom stresses are approximated by the Chezy 
relationship: 


Tv 


ci? 


( 5 . 5 ) 


where the Chezy coefficient, is given by 



iJd_(Vi + 0 
N 


( 5 . 6 ) 


where (h + ^ ) is in feet, and N is the Manning factor. 
Choice of a value for the Manning factor must be made 
experimentally. The surface stresses due to wind are 
approximated by the quadratic law for turbulent flow. 

The final differential equations used are the 


following: 
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+ ? - MM 4 - P-MV = - a 2^ 4 -f V+ kj A l WyjWy I - q U 

at ax ay ^ ax z^’ h ^ c^h 

(^•7) 


ly.+ ^ ^ _ ^--fu+iaP-. WvlWvl - qV(u 

a-b ax ^ J ay c[fH 

(5.8) 

-^-L 4- .Pjiy + ^HV =0 (s.'i) 

at 5x Dy 


where capital rj and V signify vertical averaging, as 


U=-jq- u dz av\d V=T 


■th 


H 


Vc|z 


•<-h 


is the density of air, 

kj , a dimensionless drag coefficient, is taken 
to be? 0.0025 

, fc’y are wind velocities in the indicated 
directions, 

H = h + ^ . 


The solution approach used by Leendertse involves 
space- and time-staggering of 0 and V, and a multi- 
operation computation. A concise description of this 
method, given by Hess and Whit« (18), is quoted here: 


The time step is split into two halves, and the 
time derivative taken! over the half time step. . . . 
Tn the first half time step, values of IJ and ^ are 
computed implicitly along a grid row in the 
r-direction at the time (t Then V is 

computed at the same time level explicitly. In the 
second half time step, V and ^ are computed 
implicitly at (t * 1)AT along grid rows in the 
y-direction, after which tl is calculated explicitly 
at (t + 1)AT. 

"In the first half time step, the time 
derivative of n in the x-momentum equation is 
approximated by a backward difference: ... In the 
second half time step, a forward difference is used: 
. . . Thus, over a full time step, the time 
derivative is a central difference ... 


The complete finite-difference equations (19), and 
details of the solution, are given in the appendices of 
the report by Hess and White (19). 


Application to the Providence Piver 


POCK . ..! 


The model area is to be that part of Narragansett 
Bay called the Providence River. It is actually a 
partially- mixed estuary with its circulation dominated 
by tides, although wind and gravitational circulation 
may be important at times. The fresh-water inflow 
during one tidal cycle is about six per cent of the 
tidal prism. There are throe water boundaries: the 

narrows at the mouth of the Seekonk River, at the north 
end; the mouth of the Pawtuxet River on the west shore; 
and the interface with the lower bay, which is a line 
between Conimicut and Wyatt Points. There are also two 
small rivers, the Moshassuc and the Moonasquatucket, 
which enter near the mouth of the Seekonk. The magor 
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bathymetric feature is the farty-foot ieep shipping 
channel which stretches the entire length of the moiel 
area. Figure 5.1 shows the location of the model area 
with respect to Narragansett Bay and Rhode Island Sound. 

The x-direction for the finite-difference grid 
netwoi.c is chosen to follow the east shore of the river, 
approximately the same direction as the expected mean 
flow. Since this shore is fairly straight, the matching 
of the shoreline by sguare grids is optimised. Thus the 
x-axis is a line directed S 22° E, the m-index increasing 
to the south. The y-axis is perpendicular with the 
k-index increasing to the east. The first comoutational 
field was taken as 32 by 12 grids, with a Ax of 1221 
feet. This was increased to 51 by IS (^x = Ay = 750 
feet) , in an attempt to improve the resolution. The 
most difficult geometric location to model accurately is 
the narrow mouth of the Seekonk River, which spans about 
750 feet where it merges with the model area, but is 
much narrower gust to the east. An attempt to model the 
narrowest entrance with a smaller Ax would have 
increased the storage and time requirements greatly, 
without providing useful resolution in any other part of 
the model. Figure 5.2 is a mao of the model area, 
showing the land and water boundaries. 

The depth field was orepared from O.S. Coast and 

Ceodatic Survey Chart Mo 278. The mean low-water 

soundings closest to the southeast corner of each grid 

fm k +:^) were placed in a matrix, adding the 
' z ■ 2 
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FIGURE 5.2 Providence River Showing Orientation of Grid System 





cl if f ecencp? between mean sea level and mean low water. 

An exception to this procedure was made at the Saekonic 
River boundary where the depths were chosen arbitrarily# 
to ensure that the model represented the actual 
cross-sectional area of the interface. 

The boundary condition at all land-water interfaces 
is that the velocity component normal to the boundary is 
zero. The water boundaries offer a choice of specifying 
either the tide height or the current velocity at every 
step. ®'or a river boundary with an approximately 
constant flowrate, it is easiest to specify the 
velocity, which is the flowrate divided by the 
cross-sectional area, and which will vary inversely with 
the rise and fall of the tide. The Pawtuxet River 
boundary is handled in this manner. 

Since the flowrates at the wide lower boundary are 
unknown, the tide height is specified. Good information 
is available. From tide-height histories taken by the 
Coast and Geodetic survey at Newport, Bristol, and 
Providence Harbor, Hess and White obtained the 
amplitudes and phases of the main harmonic constituents 
of the astronomical tide. The tide height as a function 
of time (1H) is: 

^ (t) = 1 -f. (*'> * (Vo (^■'°) 

n 

where n is the number of the constituent, 

f^ (t) is a function of lunar position which 
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iRodifi'^s the amplitude, 

is the amplitude, 

is the angular speed in degrees per hour, 

Vo + n is the eguilibrium argument at t=0, 

is the epoch, relative to Greenwich, England, 
t is the time in hours after reference time. 

The seventeen largest constituents are used. Their 
amplitudes at the lower boundary are obtained by 
interpolation from the known amplitudes for Newport, 
Bristol, and Providence Harbor. 

The Seekonk Hiver boundary presents problems, not 
only because it is narrow, but also because the flow is 
still tidally dominated. Since the flow through, this 
boundary reverses direction with the tides, it is 
difficult to specify the velocity. In the Narragansett 
Bay model, Hess and White handle a similar problem bv 
expressing the velocity as a function of the three lunar 
constituents of the tide. The three components of the 
flowrate are obtained by data analysis, and the total 
flowrate as a function of time is given by (18): 


< 1 = I «s (t-T,)] (5.11) 

p 

where Tw is time to first flood after high water. 
Sufficient data was not available to apply this method 
to the Seekonk boun.dary, but assuming the lunar 
influence to bo similar, the same constituents were irsed 


to calculate the flowrates there 






Specifying t.he tide height seems a more likely 
alternative, but would mask the actual river inflow. 

The most accurate method is apparently to add to the 
model grid an area representative of the area of the 
tidal Seekonk, and let it interact with the original 
model basin. This permits adding the inflow of the 
Blackstone River at a location removed from the narrow 
boundary. The computational time and storage are 
increased, but the requirement is minimized by fitting 
the extra grids as compactly as possible into extra 
space, since the actual circulation in the Seekonk is 
not of interest, it is only necessary to represent the 
storage volume accurately, and the shape approximately. 
The result is that the tidal flow and the river flow are 
both modeled satisfactorily. 

Two considerations determined the length of the 
time step. First, it must be compatible with the 
required time step of the mass transport model. The 
second consideration is Leendertse’s parameter of 
accuracy { 19) : 




(5.12) 


where h is the maximum water depth. j3^ must be of the 
order of five or less for a solution of acceptable 
accuracy. Thus, 


= fe)(750) 


\0S sec. 




=. /3 ^ X 

"‘-wax 7-". — 


A lOO-se^onrl time step turned out to ho perfectly 
compatible with the requirements of the water quality 
model . 

Verification of the Tidal Model 

There is only one set of current data that has been 
taken in the model area. Haight, in 1930, made 
raeasur'=‘mants at many stations in Narraqansett Bay, at 
all staaes of the tide (44). Six stations within the 
model area ware occupied. The magnitude and direction 
of the current, the time relative to high tide at 
Newport, and the estimated component due to wind ware 
presented . 

To compare the model predictions with tha existing 
data, runs of the model were made for a data on which 
the tidal range was the average, 4.6 feet at Providence. 
Initial conditions were established by a ?4-hour 
simulation starting with zero tide height and 
velocities. This was found to be an adequate startuo 
time to eliminate transients. 

"^he velocities calculated are dependent upon the 
bottom rouqhnass, and thus upon the value chosen for the 
Manning factor in Equation 5.6. Hess and White (18) 
recommend a value of 0,026 for the Providence River, 
based on their modeling. Tha data gathered by Haight 
presents an opportunity to adjust the Manning factor. 
Figure 5.3 shows the measured values at the narrowest 
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part of the mouth of the Seekonk, at each hour of an 
average tic^al cvcle. (The estimated components r!ae to 
wind were removed by Haight.) Using a single-grid 
boundary of the same cross-sectional area as the real 
channel at Haight’s station, the Hanning factor was 
adjusted until the best fit was found. This occurred 

for a Hanning factor of 0.030. 

some distortion is expected at this location due to 

the fact that the model grid is about three times as 
wide and ona-third as deep as the real channel. This is 
all right at mean tide, but at low tide, the 
cross-sectional area of the model is considerably lass 
than that of the real channel, due to the lesser depth. 
However, since the greatest flow occurs when the tide 
level is near mean, this effect is not severe. 

Figure 5.3 shows excellent agreement in phase and 
in amplitude with Haight’s data. The first peak of the 
double flood (characteristic of Narragansett Bay) is not 
well-matched in amplitude, but failure to reproduce a 
single data point does not indicate a flaw in the 
equations or boundary conditions. Indeed, were Haight’s 
points connected into a continuous curve, the area 
underneath would indicate that the volume of the flood 
tide is greater than that of the ebb, river inflow 
notwithstanding. The modal is checked for mass 
conservation at all times, bv summing the river and 
tidal inflows over time, and comparing this to the 
change in the water content of the entire modal. The 


greatest error found, which always returned to zero 
within the tidal cycle, was 0.7 per cent of the total 
mass. 

In an effort to verify the model more broadly, 
current vector plots were obtained for the whole field, 
at each hour of a tidal cycle. Haiqht's normalized 
vectors were transferred to the plots for comparison. . 
These are presented in Figures 5.4 through 5.15. At 
slack tide the currents are variable and agreement is 
not good. However, the model appears to reproduce the 
data eYtremely wall on the ebb tide, both in magnitude 
and in direction. Motion becomes random again at slack, 
but returns to fair, though clearly not as precise, 
agreement on the flood .tide. This situation is 
complicated by the characteristic double flood, which in 
effect inserts an extra period of slack water between 
the two parts of the flood. On the whole, it appears 
that the model predicts the same kinds of motion-- 
eddies of the same size, duration, and location-- as are 
indicated by the field measurements. The one consistent 
difference is that the measured magnitudes are greater 
than those predicted. This is very likely to be the 
result of vertical averaging, since Haight’s 
measurements were taken in the upper layers, mostly by 
floating spars. 


Input for the Water Quality Model 
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figure 5.4 


comparison of Field Measurements and Model Predictions for Tidal 
Velocities, High Water at Newport, R. I. 
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FIGURE 5.11 



comparison of Field Measurements and Model Predictions for Tidal 
VelocitieSf Seven Hours after High Water at Newport, R. I. 











figure S.14 



Comparison of Field Measurements and Model Predictions for Tidal 
V'^locities. Ten Hours after High Water at Newport, R. I. 




figure's. 15 Comparison of Field Measurements and Model Predictions for Tidal 

F.ir>vi>n Hours after High Hater at Newport, R. I. 
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ValuRs of the u- anil v-vplocities, and the tide 
heights generated by the tidal model» are stored for 
each grid point and time step. The information is read 
as input to the water quality model at each time step. 
Care must be taken to supply the information correctly, 
since the two models are computationally incompatible. 

The problem is mainly that the two-dimensional 
model has a time step divided into two levels, while the 
three- dimensional model has three levels. The 
information generated by the tidal model is u, v, and 
the tide height (called se) at each step. The 
information required by the water quality model is u and 

V at the beginning and and of each step (u, v, up, vp) , 

and the tide height at each fraction of a step and at 
the end of the step ( 1 , 1 * 5 , 1 I * 1 : se, sep, sag, 

ser) . The match-up is made by staggering. the reading of 

V and u, reading the unwanted arrays into a dummy 
variable that is not used. The order of reading and 
dummying the velocities is reversed after each step, to 
prevent the models from building errors due to model 
matching. The values read at 1 *1 are carried over to 
the values of 1 at the start of the next step. This 
sequence is illustrated in Figure 5.16. The values read 
and used are circled; the values read into the dummy arc 
crossed off. The time levels assigned to each variable 
are indicated by the subscripts. The method Imposes the 
requirement that the water quality model have a time 
step three times as long as that of the tidal model. 



Time steps, water quality model 


Figure 5.r6. Reading sequence for hydrodynamic input 

to water quality model 







Fortunately, for this body of water, it is possible to 
use the optimum time step for both moiels. 

All of this is accomplished by a subroutine in the 
water quality model which is called at each time step. 
The velocities, read into a horizontal matrix, are then 
extended to all the vertical levels. This methoi has 
been checked for mass conservation, as described in 
Chapter VIT, and found to be excellent in this regard. 


VI. SI”NT^ICANCH AND BEHAVIOR OP COLIFORH BACTERIA IN 

SEA HATER 

Valua as an Indicator of Sewage Contamination 

Due to the difficulty of isolating pathogenic 
bacteria and enteric viruses from water and sewage, it 
has long been the practice to infer the guality of 
water, or the potential hazards of wastes, from the 
concentrations of the more abundant and easily 
detectable coliform bacteria group. The groundwork for 
this practice was laid by Escherich in 1885, who 
determined Raccilus coll to be characteristic of the 
feces of warm-blooded animals. Although the total 
colifor-T! count is still the commonly used indicator and 
is the basis of water quality standards, its use as an 
indicator has recently come under heavy attack from 
microbiologists. 

Concern over the presence of disease organisms in 
natural waters falls into three areas: transmission of 

disease through drinking water (not a concern of this 
paper), contamination of shellfish, and infection of 
swimmers. Coliform standards have lono been in effect 
for the first two considerations (45,46). There is 
clear evidence linking consumption of contaminated 
shellfish or drinking water to outbreaks of typhoid 
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fever or infectious hepatitis (47,48). 

The danger of bathing in contamihiated sea water is 
not at all clear. There is almost no evidence 
conclusively linking polluted beaches to disease 
outbreaks. on this basis, some have gone so far as to 
state that bathing in sewage-polluted sea water carries 
negligible risks to health, even if the water is <^ALITY 

aesthically unpleasant (43). However, it is extremely 
difficult to detect disease contracted by swimming, lue 
to such problems as the transience of swimming 
populations, the long incubation periods (49) of 
diseases such as hepatitis (a month or more), and the 
scarcity of enteric diseases in the populations of the 
united states and the United Kingdom, where most of the 
studies have been attempted (49). This difficulty has 
been used as an argument both for and against tha 
imposition of microbial standards for beaches (49,50). 

The argument against classifying bathing waters by 
coliform levels centers upon, first, the absence of 
evidence of disease transmission, and second, the large 
observed variations in coliform counts with time at a 
given beach, which would seem to preclude assigning a 
beach to a certain class. Shuval (49) , however, 
concludes from a mathematical estimate of the 
probability of contracting disease, that standards are 
needed. Although enteric disease has not been linked to 
contaminated bathing waters recent studies point 
conclusively to the danger of skin and upper respiratory 



tract infections (51,52). It may be safe to conclude 
that when going swimming, one would wish to Icnow whether 
the water is polluted. 

The presumption that the total coliform level can 
indicate how polluted the water is, and with what 
harmful organisms, is increasingly in doubt, although 
the occurrence of pathogens such as Salmonella and 
Streptococcus is generally found to be related to the 
coliform count (47,53), this is not always the case. 
Disease outbreaks due to Salmonella and Shigella have 
occurred in instances where the drinking water met the 
coliform standard (lass than 2.2 per 100 ml.), and the 
ratio of Salmonellaa to coliforms, usually very small, 
was greater than one (52). Another problem is that 
coliform levels can increase enormously in the presence 
of organic nutrients, while this effect is not observed 
for pathogens or viruses (51,52). Furthermore, and 
perhaps the most damaging to their indicator status, 
coliforms are more susceptible to disinfection than 
enteric pathogens (52,54,55). This would mean that the 
coliform test overestimates the effectiveness of sewage 
treatment, and the quality of the receiving waters. The 
nature of the coliform group itself presents soma 
problems. Not all members of the group are of fecal 
origin; some occur on plants and in soil, and would be 
present in large numbers in runoff that was not 
necessarily contaminated with sewage. The different 
members Pscherichia coli . Klebsiella , Citrohacter, and 


Enterohacter have been found to have different die-off 
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rates and different responses to nutrients (56), The 
various methods used to obtain coliform counts do not 
even detect all of the same organisms (52). 

Dutka (52) soecifies four criteria a good indicator 
should meet. They are: 

1. Occurring in much greater numbers than the 

ORIGINAL PAGEJS 

pathogens; OF POOR QUALITY 

2. Not proliferating relative to the pathogens; 

.1, Being more persistant through disinfection and 

in the environment than pathogens; 

4. Yielding an unambiguous identification. 

Dutka concludes that the coliform group fails all these 
tests. He recommends that fecal coliforms, together 
with fecal streptococci, be used as an indicator 
instead. rt has been found that enterococci are not 
subgect to the growth phase, and the death rate is 
smaller and less sensitive to the environment than that 
of conforms (55) , Others have recommended Escherichia 
con, which is of unquestionable fecal origin and has 
been studied individually to determine its die-off rate 
(51,53). 

Despite the availability of this information, it 
will probably take years for regulatory agencies to 
adopt better indicators and acquire the new techniques. 

Meanwhile, the deterministic water quality model comes 
into its own. Any constituent can be modeled, including 
pathogens and viruses as well as indicators, providing 
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source levels and die-off rates can be estimated. Some 
good information is available on# for example# 

Salmonella (53). The use of a model can eliminate 
problems such as the uncertainty in the widely-used NPN 
(most probable number) test, and the presence of 
conforms of non-fecal origin. Field data may be used 
to verify the model under known conditions, and the 
model will describe the effects of variations in the 
conditions. 

perhaps it is fortunate that total coliform data is 
the kind most likely to be available for verification 
purposes. since much literature on coliform kinetics is 
available, this permits the best possible formulation of 
the kinetics in the model. 

Reaction Kinetics of Coliforms in Sea Water 

Tt is questionable whether the disappearance of 
coliforms in sea water is correctly called either 
"die-off" or "mortality". Inactivation and 
sedimentation are likely to be mechanisms of 
disappearance. It is clear from many studies (51,54, 
55,57,58,59) that the disappearance is much more rapid 
in sea water than in fresh water. Many studies 
attempting to estimate the decay (disappearance, 
die-off) coefficients have been made, and many differing 
results have been obtained. At this point, it is still 
uncertain which mechanisms prevail, and under what 


conditions 
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The behavior of coliforas in sea water consists of 
up to three phases. Figure 6.1 illustrates the 
different kinds of behavior that may be found. The 
first, which is not always observed, is a lag phase, in 
which the population does not decrease, and may increase 
if nutrients are present. Lag periods of the order of 
0,4 day (57> , one day (56), and three weeks (51) have 
been reported. This phase is followed by an exponential 
decay. Hany investigations have attempted to determine 
the coefficient of this decay, and the conditions upon 
which depends. Finally there is a resistant phase, 
in which a certain portion of the population will 
persist long past the rest, due to an inherent ability 
to resist the pressures of the environment (56) . This 
phase has not been well characterized. 

The exponential decay coefficient, Kj , is defined 
in terms of the ratio of the coliform count at any time 
t to the initial count (59); 



Values of are standardized by using the time, t,^^ , 
required for ninety per cent of the coliforms to 
disappear, so that 
0.1 = 

Thus, 



t .^0 


(fc.Z) 
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The quantity tcj^ is of great interest because it gives a 
guide indication of the degree of self-purification of 
which the water is capable. Values for sea water can be 
as small as 20 minutes, or nearly 200 hours (60) . 

Studies of the decay rate have yielded widely 
varying results depending upon the location, whether 
samples were stored in the lab or in situ, upon the 
counting method, and on whether artificial or natural 
sea water was used. To obtain numerical values of the 
decay coefficients, and a model of the processes, 
reviews of the past findings must be made. Fortunately, 
two recent reviews offer an opportunity to neatly 
resolve part of the coliform kinetics problem. 

ilitchell and Chamberlin (61) have formulated a 
model incorporating the major known, or generally 
accented contributors to coliform disappearance and 
growth. Disapoearance is due to sedimentation, solar 
radiation, predation, and physicochemical effects 
(osmotic effect, pH, specific ion toxicity). Growth is 
due to the presence of nutrients in the plume (carbon, 
nitrogen, phosphorus) . The mass transport equation is 
given for one-dimensional flowi 


U - V 5 dC 


d X 


dz 


=. ■+ e -t -Sin - a 

W l.Ks+5^ J 


C 


Xwt Pp 


1_ (kp + C)YpcJ 


c - - rC 


(63) 




where = sedimentation velocity; 

C = conform concentration; 

/4y^= maximum coliform grovrth rate; 

-'s ~ half-saturation constant for: 

. 3 , the concentration of nutrients; 
a = endogenous respiration rate; 

\= maximum predation rate; 

?p = concentration of predators; 
iCp = half- saturation constant for P ; 

= yield of predators on bacteria; 

=die-off rate due to solar radiation; 
l(t) = solar radiation intensity; 

«(5 = attenuation of light in water; 
r = physicochemical die-off rate. 

The trouble with this model is that most of the 
above parameters are unknown, x^ot only are there no 
general values, but there aren’t even good estimates or 
field measurements of, for example, the concentrations 
of predators. The model is presented in a more useful 
form in the following table, which gives estimated 
maximum die-off rates due to each component (6l), 
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Table 6.1. Components of the Mitchell- Chamberlin 

Kinetic Model 

wani-ftr Sensitivity K (max) 


Sedimentation 
Sunlight 

Predation 
Nutrients 

^ Physicochemical 

Still, it is not possible to make good estimates of 
the components of K^, within the given maxima, without 
extensive measurements. However, from the literature, 
it is clear that three factors dominate the behavior of 
conforms in sea water; temperature, sii^^ight, and the 
presence of nutrients. 

The temperature dependence has been observed as 
early as 1956- (57). 'fhe familiar rule of thumb for 
biochemical processes, that the rate doubled for a 
temperature Increase of 10* C, has been verified by 


degree of treatment 0.6 /hr. 

turbulence 

season /hr. 

latitude 

turbidity 

temperature 0.3 /hr. 

temperature -0.6 /hr. 

degree of treatment 
organic pollutants 
■hem-nerature 0.15 /hr. 
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aameson and Gould (60) for oollfom decay; In fact, they 
found the factor to be 1.97. -?rou a combination of 
earlier results and their own Investlsatlons in which 
light was excluded, Gamescn and Gould proposed a 
relation for t,. as a function of temperature, In the 
absence of solar radiation. 


1op:„ t,, (dark) = 2.292 - 0.0295 T 


(6.4) 


where T In in degrees Centigrade. 

She effect of sunlight, when clear beakers of sea 

,®ter were exposed to It, «s to reduce t,» to 

as 20 minutes (60). Many Investigators have observed 

that the effect of solar radiation is pronounced. 

Mdlng the effect of sunlight to Equation 6.4, then, 
should incorporate two of the three Important variables. 
Kitchen and O'namberlin obtained their estimated maximum 
Kj due to sunlight of 4 /hr. from an estimated minimum 
of 30 minutes; using Gameson and Gould's 20 minutes, 

the maximum becomes 6.9 /hr. 

In using a finite difference model, it becomes 

practical tc recalculate the decay coefficient at each 
time step, as the altitude of the sun varies. Assuming 
the maximum K of 6.9 to occur at the surface, rtth the 
sun directly overhead, the reduction In radiation 
intensity at the surface and at all depths, at other times, 
can be calculated. Halations giving radiation Intensity 
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at the surface as a function of latitude, time of day, 
time of year, surface reflection, turbidity, atmospheric 
transmissions cloudiness nave oeen set forth oy Ryan 
and Stolzenbach (62). As light is attenuated by tne 
water, now varies with depth, and a value must be 
calculated for each model grid point. Values for the 
8d:tenuation coefficient have been obtained for tne 
Providence River (63). 

Modeling the lag or growth phase is more difficult. 
The lag phase is^ usually characterised by a time period, 
on the order of a da^, during which no decay occurs*. 

This would he very difficult to incorporate into an 
Sulerian calculation, since the time a particle has 
been in the field is not known. Equation 6-3 suggests 
that nutrients could be modeled as a second constituent, 
which would make good use of the model's capabilities 
as described in Chcq>ten II- The concentration of 
nutrients • and the coefficients of interaction with 
coli forms', are: not known for the Providence River* end 
would be difficult to obtain. However, if in any 
problem^ modeling the growth phase was of particular 
interest, it might be worthwhile to attempt to obtain 
such data. This project will only attempt to 
demonstrate that the model has this capability. 




VTT. ESTUARY APPLICATION OF THE HODRL 

ORIGINAL PAGE 

OF POOR QITAii'l 

Water Quality in the Providence River 

The Providence River is surrounded by areas of high 
population density. It receives, either directly or 
indirectly, the affluent fron sewage plants serving a 
population of about 373,000 (64), plus untreated waste. 

The Providence and East Providence sewage treatment 
plants discharge directly into the model area. The 
other significant quantities of waste enter by the 
Pawtuxet and Seekon)c Rivers, 

The levels of pollution are such that the entire 
model area is permanently closed to shellfishing. Field 
measurements of the coliform levels were made in 1966 at 
about eight stations in the model area, and at the 
mouths of the Pawtuxat and the SeeJconk, by the NPN 
method. since then, the quality has been improved, 
mainly by the addition of secondary treatment at the 
Blackstone Valley sewage plant on the Seekonk River, 
the main source of sewage pollution now is the 
Providence sewage treatment plant at Fields Point. 

Since the storm and sanitary sewers in the Providence, 
Pawtucket, and Central Falls area are combined, heavy 
rainfall causes the flow to exceed the plant's capacity. 
Untreated waste is then discharged directly into the 


river. The official policy is to close nearly 10,000 ‘ 

additional acres south of the model area to 
shelltishing, for seven days after ona-half inch of 
rainfall, and for ten days after an inch or more of 
rainfall (65). The Pawtuxet River also remains a large 
source of sewage pollution. 

Figure 7, 1 shows the field measurements of total 
colifortn MPR made in 1966. Concentrations are plotted 
on a log scale, against the distance from the mouth of 
the seekonk, along an approximately central axis of the 
model, surface samples were taken at four different 
stages of the tide at each station: high tide on August 

31, low tide plus three hours on September 8, high tide 
plus three hours on September 28, and low tide on 
November 21. All samples were taken in the mornings. 

The different dates make it difficult to distinguish 
seasonal from tidal variations. However, since the Host 
Probable Number is only a statistical estimate, and not 
3 count of the coliform population, the uncertainty in 
the data itself may be greater than the seasonal 
varia tion . 

The coliform levels for the four sources to be 
modeled (Seekonk River, Pawtuxet River, Providence 

treatment plant. East Providence sewage treatment 
plant) are obtained from data kept by the Rhode Island 
Division of Water Supply and Pollution Control (64) . 
Figure 7.1 shows four counts at the mouth of the 
Seekonk. (The other, higher counts shown at L =0 are in 


H- high tide 8/31/66 
B- middle of ehh 9/28 
L- low tide 1 1/21 

P- middle of flood 9/8 


EH e 


10 STATIONS 14 17 


20 2 


L, thousands of feet: distance from mouth of Seekonk 
Figure 7.1. Field measurements of total coliform MPN 
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the moath of the very small, non-tidal Providence River, 
whose coliform contribution is neglected.) At low tide, 
ebb, and flood, the value is *»300 coli /100 ml.; at high 
tide it is 4300. Extensive discussion of the modeling 
of this boundary follows. The Pawtuxet River, being 
non-tidal, will be modeled by a constant coliform level. 
From several measurements made at different times, the 
value 4100 /100 ml. is the most often repeated, and is 
selected as the most representative., counts for the 
treated effluent of both sewage plants average 2300 /100 
ml. rising these as source levels, verification of the 
model will be attempted for the conditions of 19fi6. 

nodelino of a Conservative Constituent 

The mass transport model is now used to model the 
distribution of a constituent equivalent in source 
levels to conforms, but with no decay specified yet. 

The model grid and depth field are, of course, identical 
to those of the tidal model already developed, except 
for the additional area used by the tidal model for the 
Sechonk River. The horizontal dispersion coefficients 
are calculated from Equations 3.29 and 3.30. The value 
of is taken as a turbulent diffusion coefficient, for 
a length scale equal to a typical width of the estuary. 
From Equation 3.13, a value of 11^=30 ft* /sec. is 
obtained for 1 =3000 feet. 

A salinity field is needed to calculate the 
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vertical diffusion coefficients. The Providence River 
has been found to be highly stratified. Rather than 
doubling the computational time by modeling salt 
dispersion, a constant salinity field is obtained by 
averaging field measurements (66) over a tidal cycle, 
and estimating a linear salinity increase from the 
inshore to the open end of the model. A different 
equation is used for each level. Thus the salinity 
varies from 14 ppt. at the north end to 22.5 at the 
south end in the top level, and from 27.0 to 12.9 opt. 
in the bottom level. 

The effect of the stratification is to suppress the 
effects of turbulence and small waves. A base value, 
similar in purpose to was set at 0.001 ft /sec. 

based on the order of magnitude of the smallest measured 
values of found in the literature (.13) . As 
calculated by Equation 3.19, Oj, exceeds this value only 
when waves such as would be generated by a sustained 
20-knot wind are specified, and then only in surface 
waters. The stratification of the Providence River is 
thus seen to be a very important factor in suppressing 
vertical exchange. 

-Tt is desired to use a time step three times as 
long as that of the tidal model, or 300 seco'ids. 

Checking the velocity and dispersion criteria, the 
maximum velocity which occurs is 2.3 ft /sec, but 
velocities greater than 1.0 are rare. The largest 
dispersion coefficients obtained, prior to the use of 
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the subroutine for discontinuities, are about 80 ft 
/sec. Thus, 

At Cway) ^ 2.3 >x300 _ OPO. O r - = 0.4 

AX 750 "750 

X>v ^ aOx^OO _ 0.043 

(ax) 2 C7 50y 

A 300 -second time step is acceptable, although it 
strains the limit at the Seekonk River boundary. It i 
not deemed profitable to double the computational time 
just to increase the accuracy at this location. 

It is essential that the model conserve the mass of 
all constituents to within a few per cent, and that 
errors in mass conservation do not increase with time. 
The comouter program checks for such errors at all 
times. Another method of checking the model for mass 
conservation and stability is to model a uniform field 
of an arbitrary conservative constituent, with constant 
and equal boundary concentrations. This also indicates 
how well the method used to link the tidal and mass 
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transport molals, as dascribed in Chapter V, conserves 
mass. This has been done for a field and boundaries of 
5.0, simulating 24 hours vfith complete hydrodynamic 
input. Any deviation from the value 5.0 indicates 
computational error. The greatest deviation observed at 
any point was 0.073 (1.46 *) , and the maximum total mass 
conservation error was less than one-tenth of one per 
cent. This test indicates very satisfactory performance 
by the model. 

The ooint sources at the two sewage treatment plants 
are modeled by a source term existing in one grid or one 
column of the inolal. Due to the stratification present, 
and the fact that both outfalls are in shallow water (5 
to 10 feet) , it is expected that the effluent will rise 
to the surface. Therefore, the sources are placed in 
the surface level, N =2. The source levels are 
estimated by multiplying the coliform concentration by 
the discharge per second, and dividing by the volume of 
the grid. A flowrate of 53x10® gallons per day 
(Providence), with a coliform count of 2300 /100 ml., 
becomes a source of 0.1 ml. per second throughout the 
grid. Por the Rast Providence plant, with a discharge 
of 4.44x10® and the same concentration, the source 
strength is 0.01 /100 ml. 

Figure 7.2 shows contour maps of a simulation with 
the Providence plant as a point source. The 
constituent, with no decay specified for observation 
purposes, was injected into an empty initial field at 
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the second leval (N =3), starting at high tida. Figura 
7.2 shows the distribution after 21 hours, which was the 
middle of the flood tide. The upstream discontinuity 
problem is evident in the pronounced dip just below 
(south of) the source grid. (The flood tide is flowing 
up the page.) The variation in the vertical direction 
indicates the magnitude of diffusion for 0 =0.001 ft^ 
/sec.. The persistent 5-contour south of Fields Point, 
and the 20-contour opposite the source, are in shallow 
areas. This indicates that the coliforms become 
uniformly distributed in the column, by diffusion, when 
the death is of the order of six feet or less. It can 
also be seen that the sewage treatment plant is a small 
source of fecal contamination under normal operation. 

It appears that its contribution wiU be negligible 

compared to the river sources. 

rinder the assumption that the sewage-contaminated 
fj.pgh-water inflows will be buoyant, the source levels 
at the river boundaries are taken as maximum at the top 
level, decreasing linearly to one-third maximum at the 
bottom. Whether this distribution persists downstream 
will be an interesting facet of the three-dimensional 
model. Adding the river sources reveals several 
difficulties, due to the fact that the flow reverses at 
the mouth of the tidal Seekonk, First, specifying a 
constant boundary concentration of 9300 results in the 
upstream discontinuity problem again during flood tide. 
The concentration lust inside the boundary takes a dip. 


with the result that a large mass conservation error, ud 
to 3b par cent, is prodacaS, This is because the 
transport out of tha model is calculated with the 
concentration just inside the boundary, which is too 
low. Tt is therefore necessary to extrapolate the 
boundary concentration. The first-order extrapolation 
is given by 

_ Vk-i (Pk- Pk-i') 

^v< 

— Py ( Pk-I - Pk) C? l) 

This is used to calculate the boundary concentration 
during flood tide (whenever V at the boundary is 
positive). When tha tide reverses, the concentration 
reverts to the constant 9300 level. 

This still leaves a mass error of as much as nine 
per cent, due to the abrupt change in concentration upon 
returning to abb. This is unrealistic, because the 
Blackstone Valley sewage plant is about three miles up 
the seekonk. The flood tide would push the polluted 
waters back from the mouth, and the coliform level would 
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,t,duallv ratarn to ao tho abb prograased. It 

is becaoKaty to use a ca.p fubction to bring tha 
boandary ooaoantration baok ap- slbce bO field 
inforeatlob la ayallable to indicate ho. long this 
should take, the length of the ranp in taken to be 30 
tine stens. l.a.. the first hours of the ebb tide. 

Shen all this retlnesent l.s sade. there remains a 
mass conser.atlon error between *3.5 and -1.5 %. (» 

positiue error is an eroess of the mass in the field 
over the sum of the initial field plus the net influx; a 
negative error is the opposite.) The error is found to 
folio, the boundary velocity in phase, as shown in 
Figure 7.3. This indicate.s that the error is due to the 
degree of approximation in the extrapolation. Using 
SOPS than one inside grid to extrapolate might reduce 
the error, but was not attempted. Since the error is 
roughly proportional to the velocity, the error might 
also he due in part to the aforementioned stability 
limit. Again, the error is not severe enough to demand 
a shorter time step. The fact that the negative error, 
ths ebb tide, is only half as large as the error on 
the flood, indicates that both factors induce error- 
the extrapolation error being added to the high-velocity 
error on the flood, and having no effect on the ebb. 

Figure 7.3 also shows the boundary concentration as 
a function of time and velocity. The concentrations 
marked by the letters E, L, F, and H are the field 
ueasurements (22, at ebb, low. flood, and high tides 





113 

respectively. The large mass conservation error at the 
start is the result of the empty initial field, 

Extended runs of the model, to bring the nodal uo to 
steady-state concentrations, revealed a peculiar and 
damaging affect of the method of adjusting dispersion 
coefficients at discontinuities. Occasionally, large 
errors would appear in the concentration field, at 
apparently random locations and times. Osually these 
errors ware damped out, but in several cases the error 
continued to grow, with both positive and negative 
concentrations five or six orders of magnitude too 
large. 

!Ipon closer inspection, it was found that this error 
began at a point where the dispersion coefficients ware 
being adjusted by the subroutine, at a minor 
discontinuity. The method appears to destroy the 
stability of the overall solution technique. Although 
the mass transport computation tends to damp out the 
error, the subroutine overcompensates bv continuing to 
increase the dispersion coefficients as the error 
increases, causing the error to spread. 

The empirical coefficients, described in "hapt^'s: IV 
and figure 4.7, ware set at the rather high value of 0.5 
at this tine. Although reducing this value might have 
solved the oroblem. It was decided to bypass the 
subroutine, to remove all possibility of another such 
error. The capacity to reduce generation of negative 
concentrations was lost, but since 


the discontinuities 
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in the fieli wsrs not Important at this point in the 
computation, the adjustment nas not greatly needed. It 
is recommended that the subroutine be used whenever 
discontinuities are expected to occur, preferably with a 
smaller value of the coefficients, such as 0.2. If the 
concentration field reveals an error at any time, the 
subroutine can then be bypassed, or suitably modified to 
account tor the peculiarities of the flow system. 

continuing with the simulation of the conservative 
constituent, some concentration contour maps are shown 
to illustrate the pattern of river influx. Figures 7.4 
throuqh 7.6 Show the coliform distribution 
(conservative) IB hours from initial conditions, at low 
tide. The surface, middle, and bottom levels are shown 
in separate plots. Logarithmic contour intervals are 
used for clarity and to show the far field better. The 
most striking fact revealed here is that the vertical 
variation at the shallow mouth of the Pawtuxet River, 
from 4100 /100 ml. at the surface to 1410 at the bottom, 
is quickly eliminated by diffusion. Vertical variations 
at the deeper and more rapid Seekonk persist far 
downf ield . 

Fiinres 7.7 thcjoqh 7.9 sho« the fllstribatlon at 
23.97 hours, appto.ohinu high tiae. The poUutaht has 
bees ariuen bauh up the estuary by the flood tide. The 
uulforelty at the Pa.tuxet has disappeared. Indlcattug 
that for a diffusion coefficient of 0.001 ft /sac., 
coaplete yertlcal slxinu tales place oyer a colueu 4 






Figure 7.5. River influx, low tide 
Middle level, l8 simulated hours 




Figure 7.6. RLver influx, low tide. 

Bottom level, 18 simulated hours 



-o't 



















Figure 7.14. 
high tide. 


Steady state conservative distribution, 

Middle level, 98.17 hours, D = 0.001 ft^ 

sec. 
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Figure 7.15. Steady state conservative distribution 
high tide. Bottom level, 98.17 hours 
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feet deep, but not over the high-tide depth of feet. 
Curiouslv, the concentrations closest to the river 
increase slightly with depth. some pockets of higher 
concentration, in shallows and in coves, are left behind 
by the reversing tide, as are particularly notable in 
the lower level. 

Simulation continues, now adding the two sewage 
plant sources, which are quickly swallowed up by the 
much greater seakonk influx. After 9B hours, nearly H 
tidal cycles, material is distributed throughout the 
model area, which appears to be approaching steady 
state. Figures 7.10 through 7.12 show the low tide 
conditions at 92.25 hours, aid Figures 7.1.1 through 7,15 
show the late flood tide at 98,17 hours. It is evident 
that the vertical variations are much reduced, and are 
scarcely revealed at all by the chosen contour 
intervals, Fiaure 7.16 shows a plot of concentration 
and tido height versus time, and their clearly inverse 
relationship. As the tide ebbs, the more polluted waters 
upstream are swept past the point, increasing the 
coliforra concentration. As the tide reverses and 
cleaner water is swept back upstream, the concentration 
drops. Bullock Cove, a very shallow, narrow-necked 
inlet on the lower right, has finally filled with 
material, particularly after the flood tide. 

The vertical variation that remains is greatest 
close to the river boundaries. To show this more 
clearly. Figures 7,13 and 7.19 show concentration-depth 
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profilr*?; on linas extending away from the two river 
mouths, as shown by Figure 7,17. The rapid mixing at 
the Pawtuxet is obvious, particularly in the low tide 
case f^iqure 7.18). fluch of the vertical variation at 
the ser-konk is eliminated by passage through the varv 
deep channel, which spreads the material vertically. 

The pr'senca of the providence sewage source, which is 
located in the last grid plotted for the Seekonk influx, 
is just visible as a slight increase in concentration at 
the surface. 

Since the base value of = 0.001 ft /sec. may be 
unrealistically high, the conservative constituent run 
was reneated for a base value of 0^ = 0.00001 ft /sec. 

The profiles are shown again in Figures 7.20 and 7.21, 
revealing that the vertical structure persists much 
longer. Knowing which value is more nearlv correct 
would require field information on the vertical 
variation of coliform densities. 

ORIGINAL PAGE M 
OF POOR QUALITY 

Modeling of Coliform necay 

In order to observe the capabilities of the reaction 
matrix method, a uniform field of 50 coli /100 ml* 
throughout the Providence River was simulated with all 
velocities set to zero. The computation of depth- and 
time-varving decay coefficients, depending on the 
intensity of solar radiation, was programmed. The decay 
of the uniform field was allowed to proceed, with no 
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Figure 7.19. Concentration vs. depth near river sources. 
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ad(3itir)n of material, beginning at 6:00 A.M. on the 
longest day of the year. 

Figure 1,22 shows the concentrations at the surface 
and bottom, after 12 hours of decay, with the depth 
field nresentad for comparison. Tho decay is seen to be 
rapid In the shallow water, and apparently no more rapid 
than decay in the dark in the deepest water. It is 
evident that this rate of decay will quickly eliminate 
the coliform population in shallow water, on the order 
of 90 ^ in a day, although they may persist in deep 
water . 

The same approach was used to determine what sort of 
a reaction matrix would effectively simulate a lag or 
growth phase, taking a nutrient of undetermined nature 
as a second constituent. Figure i,22 shows the type of 
behavior produced, and the reaction matrices used. The 
values have the following moaning: 

K|, - natural decay of conforms in the dark 
K,- - growth of coliforms due to the 

(CEIGINAL PAC-l 

presence of a nutrient OF POOK QUALl* 

- d isappea rence of the nutrient due to use 
by the coliforms 

decay of the nutrients, taken as 
zero because it is unknown. 

The first simulation shows a slight lag phase and a 
rapid decav of the nutrients; the second reaction 
matrix, with a larger and a smaller , shows a 

growth ohase. It can be seen that as the population of 






coliforvns- dronr, off, ttio ra 


of non suinpf ion of tb- 


nufrisnt ieoroasas. 

Th--' Doint of this damonst ration is to show that th? 
reaction matrix can be used to model more complex 
behavior than simple decay. Some careful field or 
laboratory raeasuremant s would be needed to obtain the 
proper reaction matrix for a problem in which a lag or 
growth is expected. It should be kept in mind that the 
coliform and nutrient concentrations will probably be 
expressed in different terms, and the reaction matrix 
values must account for this. For example, if the 
arbitrary 50 -concentration used for the initial level of 
the nutrient represented 50 % sewage effluent present, 
the coefficient for coliform growth, K,^, would probably 
be much higher. A simple simulation, such as that just 
performed, would he very helpful to check the reaction 
matrix before applying it to a real-world problem. 

Figures 7.24 and 7,25 show the comparison of the 
field data with the predicted concentrations of the 
conservative constituent at the same stations (the 
surfac * values) . They appear to agree remarkably well. 
This could mean either that there is no decay or that it 
is only of the order of the difference between the 
near-steadv state prevailing here and true steady state. 
Clearly, the field data indicate no such cate of decay 
as would be produced by solar radiation. For a final 
run, a decay coefficient equal to the rate for decay in 
the dark at a water temperature of 75®F (Kguation f>.4) 












was spGcified. 


This was run start-.inq from tha 
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consprvative concentration field at the end of 98.17 
hours, for another 18 hours. Figures 7.26 ani 7.27 show 
the daha and predictions after 12.0 and 18.0 additional 
hours. I’he predicted values have now dropped 
consistently below the field measurements. 

Since the conservative values agree so well with the 
data, it is most likely that little decay actually takes 
place under the conditions prevailing, except in th« 
southernmost part of the model area, most distant from 
the sources. This appears to indicate that nutrients 
are present throughout the area, causing the coliforms 
to persist. In 1966, the Blackstone Valley sewage plant 
on the seekonk River was a very large source of organic 
matter, and the Pawtuxet River is thought to be a source 
of a variety of nutrients (56). This is the likely 
explanation for the apparent lack of decay. 

It can be seen that in order to model coliforms 
faithfully, a good idea of the processes prevailing is 
needed. This can only be obtained from field 
measurements in the area to be modeled, or in a very 


similar area 
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Figure 7.R6. Model predictions and field data, 
K.* 0.078/hr, low tide, 98.17 + I 8 hours 
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STEADY STATP! MODEL 


Since there ace many situations when the steady 
state solution of the mass transport equation is of 
particular interest in an area under study, it would be 
desirable to determine necessary modifications to the 
existing numerical procedure to handle this situation. 
Following the work of Douglas (17) and Wachpress (67), 
the time step increment At in Equations 2.20, 2.21, and 
2.22 can be replaced by a positive number iteration 
parameter or sequence of iteration parameters and by 
iteration with these parameters a steady state solution 

obtained. 

After convergence of the solution is assumed, the 
problem then becomes the determination of a sequence of 
iteration parameters which, when applied in some cyclic 
pattern, will cause the rate of convergence to be 
maximized. Since the literature (17,25,67) provides 
only an indication of possible iteration parameters for 
a simple heat diffusion problem with constant dispersion 
coefficients, an optimum sequence of parameters is not 
available for the general mass transport equation and 
normally has to be determined through numerical 
experiments. Indications of possible parameter 
selection have been male in the work of A7i7 and Heliums 
(4.1), cordon and Spaulding (68), and Anasoulis and 
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McDonald (f>9) dut ar? not directly applicable to this 
case. 

Sciiilding (24) , in earlier work on the three- 
dimensional convective dispersion model, has used a 
cyclic iteration sequence given by 

25 X 10 

1.25 X 10“ 

5.25 X lO'"^ 

3.125 X 10"^ 

for relatively simple channel flow cases. There appears 
no reason to consider this an optimum sequence, but the 
solution converges quickly with reasonable error levels, 
so this sequence will be used in the work that follows. 

nsinq the iteration sequence outlined above, th^* 
steady state model has been compared to appropriate 
analytical solutions for conservative waste discharges 
in open channels with and without decay as wall as the 
simple dissolved oxygen-biochemical oxygen demand 
coupled reaction mechanism. With an iteration 
convergence criteria of 5 x 10""^ , the maximum error 
level in any solution was of the order of 0.5-1 'R, 

To further tast the steady state model, a comparison 
between an analytical solution for a continous point 
release In a three-dimensional uniform channel flow and 
the numerical solution ware compared. Table 3.1 givas 
the details on tha model parameters employed. 

Fiquca 8.1 shows a comparison between the analytical 
solution of rlaary and Adrian (70) for concentration and 
the present model alonq a line passing through the point 

ORIGINAL PAGE IS 
OF POOR QUAUTSr 


145 


TABLE 8.1 

COOTINUOUS POINT SOURCE RELEASE IN A UNIFORM CHANNEL PLOW 


MODEL PARAMETERS 


FLOW 

DISPERSION 

SOURCE 

GRID SPACING 
CHANNEL GEOMETRY 

ITERATION CONVERGENCE CRITERIA 


UNIFORM CHANNEL 
FI,OW U = 1 FT/SEC 

D = D =500 FT^/SEC 

D^ = ,8l FT^/SEC 
z 

3 

.03 UNITS /FT SEC 

AT X = 5000 FT, Y = 4500 FT 

AND Z = 12.5 FT 

Ax = 1000 FT, Ay = 1000 ft 
Ati=Az = 5FT 

CHANNEL DIMENSIONS, 

LENGTH - (DIRECTION OF FLOW) 

30,000 FT 

WIDTH - 10,000 FT 

DEPTH - 25 FT 

-4 

5 X 10 


BOUNDARY CONDITION 


UPSTREAM - FIXED AT ZERO 
DOWNSTREAM - EXTRAPOLATED 
FROM INSIDE THE FIELD 
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sourcp in the -iownstreain lirection. Figure B.2 shows 
coirparisons for two stations on the outside edge of the 
channel it the same Z height as the point source. It is 
to be noted that the concentrations on the left-hand 
side of the channel (facing downstream) are higher than 
thos=> on the right. This is a result of the 
3 ff_co!itered position of the point source (y = 4500 feet 

instead of y = 5000 feet) . 

Figures 8.1 and 8.2 clearly show that the steady 
state numerical model accurately predicts the 
concentration distribution for the case under 
consideration. The one area where model solutions 
deviate significantly from the analytical solution is 
near the point source. This difference is easily 
understood since the grid spacing of the numerical 
scheme is not refined enough to represent the steep 
concentration gradients near the point source. 

Verification of the steady state model for these 
analytical solutions gives preliminary indication of the 
numerical behavior and the validity of the computational 
scheme. It remains, however, to extend the model to 
cases for which analytical solutions are not available. 

In this light, the three-dimensional mass transport 
modal was coupled to a steady state river hydrodynamics 
model to predict the motion of pollutants injected into 
separate streams for the case of river confluence. 
Employing Leendertse's (19) two-dimensional vertically- 
averaged hydrodynamics model with the input as specified 
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FIGURE 8.2. 


Comparison of numerical model prediction and analytic solution for a steady state 
point release in a uniform channel flow, Z = 12.5 Ft., y = 500 Ft., and 
y = 9500 Ft. 
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in Tablff 8.2 and ths dapth field as shown in figure 8.3, 
the computational scheme was run until the flowrates for 
all river cross-sections achieved a steady state. 

Figure 8. a shows the resulting steady-state velocity 
vector plot for the vertically-averaged circulation in 
the river confluence. 

Employing the input conditions specified in Table 
8.2, the mass transport model was run for a point source 
release in each stream. The results of the numerical 
solution are shown in Figures 8.5 through 8.9 for each 
of thf' nondimensiona 1 levels. (Level 2 at the river 
stream bottom through level 8 at the stream surface with 
a source inout at level 4.) A well-defined vertical 
distr ibu+-ion of concentration is readily noted near the 
source input points. As one preceeds downstream the 
mixing of the two streams can be seen over the PAGE 

as well as the vertical direction. QUALITYi 

Figures 8.10 through 8.14 show the exact same case 
previously described but incorporate a first-order decay 
process with a decay coefficient of 0.00001 sec. The 
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reduction of the concentrations at all levels is noted, 
but. the structure of the concentration distribution 

remains fundamentally the same. 

Tn each of the simulation cases described above, the 
outflow boundaries were allowed to seek a level 
appropriate with the internal solution by use of a 
simple continuative boundary extrapolation approach. 

From a numerical viewpoint this adjustment of boundaries 
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MODEL INPUTS FOR STEADY STATE RIVER CONFLUENCE 


RIVER MODEL 

3 3 . 

BOUNDARY CONDITIONS - CONSTANT INPUT FLOW RATES OF 100 FT /SEC AND 200 FT / 
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(U) = 0) 

- (WHEN EMPLOYED) FIRST ORDER DECAY WITH K = .00001 SEC 
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CONCENTRATION OF CONSTITUTENT (MG/L) FOR LEVEL - 2 
IN THE X AND Y PLANE 
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FIGURE 8.5. Concentration of pollutant ( Mg /1) for level -2 
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CONCENTRATION OF CONSTITUENT ( MG / L ) - 1 FOR Z LEVEL - 4 
IN THE X AND Y PLANE 
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CONCENTRATION OF CONSTITUENT ( MG / L ) - 1 FOR Z LEVEL - 5 
IN THE X AND Y PLANE 



1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

1 

0.0 

0.0 

o 

• 

o 

o 

• 

o 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0 .0 

2 

0.0 

0.000 

0.000 

0.000 

0.0 

0.0 

0.0 

0.0 

0.000 

0.000 

0.000 

0.0 

3 

0.0 

0.000 

0.000 

0.000 

0.0 

0.0 

0.0 

0.000 

0.000 

0.000 

0.000 

0.0 

4 

0.0 

0.000 

0.000 

0.000 

0.0 

0.0 

0.0 

0.000 

0.000 

0.000 

0.000 

0.0 

5 

0.0 

0.0 

0.000 

0.000 

0.0 

0.0 

0.0 

0.0 

0.003 

0.005 

0.002 

0.0 

6 

0.0 

0.0 

0.000 

0.000 

0.0 

0.0 

0.0 

0.0 

0.003 

0.020 

0.0 

0.0 

7 

00 . 

0.0 

0.000 

0.000 

0.000 

0.0 

0.0 

0.0 

0.10 

0.06 

0.0 

0.0 

S 

0.0 

0.0 

0.000 

0.000 

0.000 

0.0 

0.0 

0.0 

0.014 

0.240 

0.0 

0.0 

9 

0.0 

0.0 

0.0 

0.000 

0.000 

0.0 

0.0 

o 

• 

o 

0.193 

0.835 

0.846 

0.0 

10 

0.0 

0.0 

0.0 

0.000 

0.000 

0.0 

0.0 

0.0 

0.0 

0.847 

0.843 

0.0 

11 

0.0 

0.0 

0.0 

0.001 

0.004 

0.041 

0.094 

0.0 

0.0 

0.843 

0.843 

0.0 

12 

0.0 

0.0 

0.0 

0.003 

0.024 

0.109 

0.020 

0.0 

0.843 

0.843 

0.843 

0.0 

13 

0.0 

0.0 

0.0 

0.0 

0.232 

1.495 

0.837 

0.0 

0.843 

0.843 

0.0 

0.0 

14 

0.0 

0.0 

0.0 

0.0 

0.282 

1.006 

0.890 

0.856 

0.843 

0.843 

0.0 

0.0 

15 

0.0 

0.0 

0.0 

0.0 

0.0 

0.991 

0.915 

0.861 

0.844 

0.0 

0.0 

0.0 

16 

0.0 

0.0 

0.0 

0.0 

0.0 

0.969 

0.920 

0.875 

0.853 

0.0 

0.0 

0.0 

17 

0.0 

0.0 

0.0 

0.0 

0.0 

0.971 

0.928 

0.880 

0.858 

0.0 

0.0 

0.0 

18 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.922 

0.893 

0.976 

0.0 

0.0 

0.0 

19 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.906 

0.885 

0.872 

0.0 

0.0 

0.0 

20 

0.0 

0.0 

0.0 

0.0 

0.0 

0.900 

0.900 

0.888 

0.0 

0.0 

0.0 

0.0 

21 

0.0 

0.0 

0.0 

0.0 

0.0 

0.898 

0.897 

0.890 

0.0 

0.0 

0.0 

0.0 

22 

0.0 

0.0 

o 

• 

o 

0.0 

0.898 

0.897 

0.896 

0.891 

0.0 

0.0 

0.0 

0.0 

23 

0.0 

0.0 

0.0 

0.0 

0.898 

0.897 

0.895 

0.892 

0.892 

0.0 

0.0 

0.0 

24 

0.0 

0.0 

0.0 

0.0 

0.898 

0.897 

0.895 

0.893 

0.893 

0.0 

0.0 

0.0 

25 

0.0 

0.0 

0.0 

0.0 

0.898 

0.897 

0.895 

0.894 

0.893 

0.0 

0.0 

0.0 

26 

0.0 

0.0 

o 

• 

o 

0.0 

0.0 

0.898 

0.896 

0.894 

0.894 

0.0 

0.0 

0.0 

27 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.896 

0.895 

0.894 

O 

« 

o 

0.0 

0.0 

28 

0.0 

0.0 

o 

• 

o 

o 

• 

o 

o 

• 

o 

0.0 

0.897 

0.896 

0.895 

0.0 

0.0 

0.0 

29 

0.0 

0.0 

o 

• 

o 

0.0 

0.0 

0.0 

0.897 

0.896 

0.895 

0.0 

0.0 

0.0 

30 

0.0 

0.0 

0.0 

0.0 

0.0 

o 

• 

O 

0.898 

0.897 

0.896 

0.0 

0.0 

0.0 

31 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.898 

0.897 

0.896 

0.0 

0.0 

0.0 

32 

0.0 

0.0 

0.0 

0.0 

0.0 

O 

• 

O 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

FIGURE 8 . 8 . 

Cone entration 

of pollutant ( Mg / 1 ) for 

level - 

5 






156 




CONCENTRATION OF CONSTITUENT (MG/L) FOR Z LEVEL - 6 
IN THE X AND Y ELANE 
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CONCENTRATION OF CONSTITUENT ( MG / L ) FOR Z LEVEL - 2 
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IN THE X AND Y PLANE , DECAY COEFFICIENT K — .00001 ( 1 / SEC ) 
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FIGURE 8 . 11 . Concentration of pollutant ( Mg / 1 ) for level - 3 , decay coefficient K = .00001 sec - 1 
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CONCENTRATION OF CONSTITUTENT ( MG / L ) FOR Z EEVEL - 5 
IN THE X AND Y PLANE , DECAY COEFFICIENT K = .00001 ( 1 / SEC ) 


0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 


0.0 

0.00000 

0,00000 

0.00000 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0,0 

0.0 

0.0 

0.0 

0.0 

0,0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 


0.0 

0.00000 

0.00000 

0.00000 

0.00000 

0.00000 

0.00000 

0.00000 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 


0.0 

0.00000 

0.00000 

0.00000 

0.00000 

0.00000 

0.00000 

0.00000 

0,00001 

0.00002 

0.00121 

0.00362 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 


0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.00000 

0.00001 

0.00002 

0.00041 

0.00417 

0.02560 

0.23077 

0.27855 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.85523 

0.85388 

0.85284 

0.85246 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 


0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.04107 

0.10779 

1.48998 

0.99806 

0.97899 

0.95298 

0.95091 

0.0 

0.0 

0.86542 

0.86227 

0.85828 

0.85563 

0.85366 

0.86240 

0.85208 

0.0 

0.0 

0.0 

0.0 

0,0 

0.0 


0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.09315 

0.20266 

0.83090 

0.87963 

0.90103 

0.90291 

0.90827 

0.89887 

0.87951 

0.87011 

0.86472 

0.86090 

0.85760 

0.85486 

0.85296 

0.85193 

0.85136 

0.85079 

0.85043 

0.85039 

0.85039 

0.0 


0.0 

0.0 

0.00009 

0.00024 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.84141 

0.84548 

0.85768 

0.85990 

0.86977 

0.85919 

0.85989 

0.85937 

0.85836 

0.85607 

0.85403 

0.85228 

0.85111 

0.95044 

0.84999 

0.84966 

0.84970 

0.84970 

0.0 


0.0 

0.00000 

0.00001 

0.00042 

0.00308 

0.00347 

0.01063 

0.01436 

0.19188 

0.0 

0.0 

0.83343 

0.83066 

0.83041 

0.82905 

0.83639 

0.83778 

0.85259 

0.84717 

0.0 

0.0 

0.0 

0,85205 

0.85124 

0.85029 

0.84957 

0.84916 

0.84893 

0.84883 

0.84894 

0.84894 

0.0 


0.0 

0.00000 

0.00007 

0.00039 

0.00572 

0,02019 

0.06854 

0.23973 

0.83413 

0.84361 

0.83851 

0.83443 

0.83240 

0.83067 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 


0.0 

0.00000 

0.00004 

0.00084 

0.00209 

0.0 

0.0 

0.0 

0.84175 

0.83699 

0.83586 

0.83328 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 


0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 


figure 8 .^ 13 . 


concentration of pollutant (Mg/1) for level - 5, decay coefficient K .00001 sec - 1 


161 














1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 


0,0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0,0 
0.0 
0.0 
0.0 
0.0 


0.0 

0.00000 

0.00000 

0.00000 

0.0 
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FIGURE 8.14. concentration of pollntant (M,/l) for level - 6, decay coefficient 
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to interior solution is extremely difficult to obtain in' 
a steady state model. The primary difficulty is that 
the unsteady boundaries tend to cause the solution to 
become non-convergen t. To overcome this, problem, the 
present study held the boundaries fixed until the the 
numerical predictions had nearly converged, and then the 
boundary was adjusted. This process was repeated until 
neither the boundary nor the internal mass distribution 
changed significantly, i.e. within the convergence 
criteria . 

Although there exists no data to compare the river 
confluence situations, the results indicate that the 
model is capable of giving reasonable quantitative 
predictions for complicated river geometries, flow 
conditions, and a variety of point loading situations. 
Extension of the model to a variety of species 
interactions is simply obtained through the reaction 
matrix- source-sink vector approach as previously 
outlined and demonstrated. 



IX. APPLICATION TO BLOCK ISLAND SOUND 

As a final application, the three-dimensional mass 
transport modal has bsan coupled to a t wo— d imansional 
vertici lly-a vera ^ad tidal model for tha Block Island 
sound area (figure 9.1) and employed to simulate a 
continuous point relaase. The emphasis in performing 
this work is to show tha feasibility of applying the 
computational system to a realistic coastal ■zone area. 

pmnloying the depth contours shown in Figure 9.2 and 
the information contained in the first half of Table 
9.1, a two-dimensional vertically-averaged tidal modal 
amploying the method of Leendartsa (19) has been 
developed for tha study area. Figures 9..1 through 9,15 
show the model predictions for zero hours through twelve 
hours after high water at Newport, Rhode Island as well 
as the tidal heights for the area (note insert at the 
top of each figure). These predictions have bean 
compared to existing National Ocean Survey tidal charts 
(71) , and the ganaral flow directions show good 
gualit.itive agreament. 

To obtain a preliminary indication of the transport 
of material in tha study area, a simple continuous point 
release of waste was simulated by a two-d imansional 




vertically-average^l concentration model (72) from a 
proposed discharge site in Charlestown, Rhode Island. 
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TABLE 9.1 


TIDAL HYDRODYNAMIC AND MASS TRANSPORT MODEL INPUTS FOR BLOCK ISLAND SOUND 


HYDRODYNAMICS 


GRIB SPACING 


TIME STEP 


A X = Ay = 6076 FT 

At = 124.2 SEC 


MANNING ROUGHNESS 


BOUNDARY CONDITIONS 


INPUT TIDAL HEIGHTS ON OPEN 
BOUNDARIES 


DEPTH 


SEE FIGURE 9.2 


MASS TRANSPORT 


GRID SPACING 


TIME STEP 


dispersion coefficient 


POINT SOURCE 


VELOCITY 


X = Ay ® 6076 FT 

/\ t = 496.8 (TWO DIMENSIONAL 

VERTICALLY AVERAGED MODEL) 

A t = 1863.0 (THREE DIMENSIONAL 
MODEL) 

D = D„ = 35 FT^/SEC 
X y 

Dg = .002 - .01 (PARABOLIC 

PROFILE) FT.2/SEC. 

4,685,776 LBS/DAY AT MID DEPTH 
(THREE DIMENSIONAL) 

FROM HYDRODYNAMICS MODEL WITH 
CYCLE USE OF THE PREDICTED 
TIDAL CYCLE 
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AFTER HIGH WATER AT NEWPORT* RI 

FIGURE 9.4 Tidal Oirrents for Block Islarid Sound in Knots One Hour After High Water 
at Newport, RI 
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FIGURE 9.5 


TIDAL CURRENTS IN KNOTS TWO HOURS 
AFTER HIGH WATER AT NEWPORT > RI 

Tidal currents for Block Island Sound in Knots Two Hours After High Water 
at Newport, RI. 
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Tidal Currents for Block Island Sound in Knots Three Hours After High 
Water at Nev^port, RI 
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FIGURE 9.6 



TIDRL CURRENTS IN KNOTS FOUR HOURS 
AFTER HIGH WATER AT NEWPORT* RI 


i FIGURE 9.7 Tidal Currents for Block Island Sound in Knots Four Hours After High Water 

at Newport, RI 
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TIDAL CURRENTS IN KNOTS SEVEN HOURS 
AFTER HIGH WATER AT NEWPORT » RI 

FIGUEiE 9.10 Tidal Currents for Block Island Sound Seven Hours After High Water 
at Newport, RI 
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TIDflU CURRENTS IN KNOTS EIGHT HOURS 
RFTER HIGH WRTER RT NEWPORT » RI 

PIGOSE 9.U Tl4al currents fur Block Island Sound in Knots Eight Hours After High 
Waiier at Newpoi^^r 












TIDAL CURRENTS IN KNOTS TEN HOURS 
AFTER HIGH WATER AT NEWPORT » RI 


FIGURE 9.13 Tidal Currents for Block Island Sound in Knots Ten Hours After High Water 
at Newport, RI 
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TIDAL CURRENTS IN KNOTS TWELVE HOURS 
AFTER HIGH WATER AT NEWPORT » RI 


FIGURE 9.15 Tidal Currents for Block Islcuid Sound in Knots Twelve Hours After High 
Water at Newport, RI 
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Figures '4.16 through 9.30 show model predictions for one 
through thirteen hours after high water at Newport in 
hourly increments. Figure 9.29 shows the concentrations 
at 5 tidal cycles after discharge began at approximately 
slack water. As expected from observing the velocity 
field, the waste cloud displays a predominantly 
along-the-shore motion that responds directly to the 
flooding and ebbing of the current in the area. The 
increase in concentration in the outfall area is easily 

seen . 

Using essentially the same input as for the 
two-dimensional vertically-averaged concentration model 
(see Table 9.1) the three-dimensional mass transport 
model was used to simulate the same release, but the 
waste was discharged at mid-depth. Figures 9.30 through 
9.34 present concentration contours at nondimsnsional 
levels 2 (at sea bottom) through 6 (sea surface) for one 
hour after high water at Newport, The vertical 
distribution of waste is clearly seen. Figures 9.35 
through 9.39 show similar plots for six hours after 
discharge began, while Figures 9.40 through 9.44 give 
the concentration distribution at twelve hours after 
waste release. The vertical stratification and increase 
in concentration with time are readily apparent, 
comparison of the three-dimensional predictions to the 
vertic^^lly-averaged case show similar behavior with 
predominant along-the-shore pollutant transport. 

Results of these three-dimensional predictions indicate 
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CONCENTRATIONS CMG/L)^ 10.0 HOUR TWO 

Concentration (mgA) Contours for Continuous Release 
by vertically Averaged Concentration Model Two Hours After High 

Water at Newport, R.l. 
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FIGURE 9.19 Concentration (mg/1) Contours for Continuous Release Predicted 

by Vertically Averaged Concentration Model Four Hours After High 
Water at Newport, R.l. 
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FIGURE 9.20 Concentration (mg/1) Contours for Continuous Release Predicted 

by Vertically Averaged Concentration Model Five Hours After High 
Water at Newport, R.l, 
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concentration (mg/1) Contours for Conti^u^s Release 
by vertically Averaged Concentration Model Sxx Hours After High 

Water at Newport, R.I. 
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Concentration (mg/1) Contours for Continuous Release Predicted 
by Vertically Averaged Concentration Model Eight Hours After High 
Water at Newport, R.l, 
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Distance (Nautical Miles) 

CONCENTRRTIONS <MG/L)f 10.0 HOUR ELEVEN 

FIGURE 9.26 Concentration (mg/1) Contours for Continuous Release Predicted 

by VerticsLlly Averaged Concentration Model Eleven Hours After High 
Water at Newport, R.I. 
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FIGURE 9.27 Concentration (rog/lY Contours for Continuous Release Predicted 

by Vertically Averaged Concentration Model Twelve Hours After High 
Water at Newport, R.I. 
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FIGURE- 9.28 Concentration (mg/l) Contours for Continuous Release Predicted 
by Verticeilly Averaged Concentration Model Thirteen Hours After 
High Water at Newport, R.I. 
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FIGURE 9.34 Concentration (mg /I) Contours for Level >.6 Continuous Release 

Predicted by the Three Dimensional Concentration Model One Hour 
After High Water at Newport, R.I. 
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Distance (Nautical Miles) 

CONCENTRRTIONS <M6/L)4 10. 0 HOUR SIX LEVEL 3 

FIGURE 9.36 Concentration (mg/1) Contours for Level-3 for Continuous Release 
Predicted by the Three Dimensional Concentration Model Six Hours 
After High Water at N/iwport, R.l. 
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concentration (mgA) Contours for Level -4 for Continuous Release 
Predicted by the Three Dimensional Concentration Model Six Hours 
After Hiqh Water at Newport, R.I. 
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FIGURE 9.39 Concentration (mg/1) Contours for Level-6 for Continuous Release 
Predicted by the Three Dimensional Concentration Model Six Hours 
After High Water at Newport, R.I, 
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Concentration (mg/1) Contours for Level-2 for Continuous Release 
Predicted by the Three Dimensional Concentration Model Twelve Hours 
After High Water at Newport, R,I, 
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FIGURE 9.41 


Concentration (mg/1) Contours for Level-3 for Continuous Release 
Predicted by the Three Dimensional Concentration Model Twelve Hours 
After High Water at Newport, R.I. 
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Sharp concentration gradients along one side of contour 
plots. This is a direct result of the point-loading 
condition and its associated depression of the 
concentration field upstream of the discharge. This 
condition could be removed by introducing artificial 
dispersion upstream of the release, but it was felt that 
the worst case should be indicated to illustrate the 
effects of the problem for realistic coastal zone 
simulations. 

Figures 9.45 through 9.49 display the threap 
dimensional results after five tidal cycles at 
approximately slack water. zomparison of these figures 
with those for the two-dimensional vertically-averaged 
case show almost exact comparison. This close agreement 
can be explained simply by noting that the 
three-dimensional model employed a vertical diffusion 
coefficient typical of a well-mixed area, and therefore 
the three-dimensional predictions should approach the 
two-dimensional case as time proceeds. The same affect 
has been noted for the Astuary application, even though 
the diffusion coefficient indicated stratification, 
because the water was relatively shallow. 

Additional complexity such as multiple-point 
time-varying loadings, multi-stage reaction mechanics, 
and stratification can be readily incorporated in the 
modeling scheme. «ith this capability it is felt that 
the modeling approach presented allows a more realistic 
prediction technique for oollutant transport in coastal 
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CONCENTRRTIONS CMG/D+ 10.0 AFTER 5 TIDE CYCLES L=2 

FIGURE 9 45 Concentration (mgA) Contours for Level-2 for Continuous Release 

Predicted by the Three Dimensional Concentration Model Five Tidal Cycles 
After High Water at Newport, R.I. 
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Concentration (tag /I) Contours for Level -3 for Continuous Release 
Predicted by th^. Three Dimension€0. Concentration Model Five Tideil 
Cycles After High Water at Newport, R.I. 
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The computational accuracy and usefulness of the 
three-dimensional mass transport model have been 
demonstrated sucessfully by the investigations and 
applications herein. The capacity to define vertical as 
well as lateral variations has been seen to be valuable, 
particularly in short-term phenomena. Over longer 
periods of time, vertical variations are much reduced by 
diffusion, at least when using a vertically averaged 
two-dimensional modal to obtain the hydrodynamic input. 

The analytical investigations of dispersive and 
dissipative effects, performed for the one-dimensional 
mass transport equation with constant dispersion 
coefficients, have been verified numerically for the 
three-dimensional model. 

The generation of computational discontinuities at 

extreme concentration gradients can be reduced, but not 

eliminated, by the method of adding artificial 

dispersion at appropriate times and places. Care must 

be taken that the increase in dispersion is not large 

enough to cause severe distortion of the concentration 

field in the area under study. The values used here for 

A 

the overall dispersion coefficients appear from the 
verification attempt to represent very well the 
processes taking place in the the Providence River. 

Two qualities of the model are particularly useful 
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in modeling of coliform bacteria distribution. The 
vertical resolution permits the incor. ’ration of the 
effect of solar radiation, and the reaction matrix can 
accomodate a lag or growth phase. As the behavior of 
conforms can vary widely depending upon the prevailing 
conditions, it is recommended that either field studies 
be made first in the area to be modeled, or the results 
of earlier studies be consulted. 

If difficult boundary conditions, such as those at 
the seaV-.onlt Piver, are seen to produce a significant 
mass conservation error, effort should be made toward 
reducing it. Elaboration upon the extrapolation 
technique is expected to help, A conclusion of this 
boundary condition research is that the maqor portion of 
the error in simulating a realistic pollutant transport 
problem is usually the error of the boundary 
approximation, and therefore more detailed understanding 
of those approximations is necessary. 

The development of a steady-state mode of operation 
has provided an additional tool to predict water quality 
for complex constant flow areas. Comparison of model 
predictions to analytic solutions and a river confluence 
case have shown the present model scheme to be accurate 
and reasonably efficient. Further testing of this 
modeling approach should be performed by comparing 
predictions for multi-stage reacting constituents to 
data for a realistic steady-state coastal zone 


circulation system 
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Coupling of the two-t^imensional vertically averaged 
tidal hydrodynamics model to the three-dimensional water 
quality model for a typical open coastal area such as 
Block Island Sound has shown the feasibility and 
capability of the present approach. Application and 
verification of this approach should build more 
confidence and experience in the ultimate use of these 
models in coastal zone management. It is therefore 
recommended that the model be verified for a number of 
coastal zone pollutant transport problems. 

In short, the three-dimensional mass transport 
model, as now developed, appears to be capable of 
handling almost any water quality problem currently 
under investigation. The major limitation is, still in 
obtaining adequate data to verify the models, and in the 
refinement of various circulation models which may be 
used to obtain the hydrodynamic input. New models are 
under development which will enhance the usefulness 
of this modal. Leendertse et al. (73) have 
developed a three-dimensional model incorporating the 
salt equation. This model has the additional capability 
of computing w, the vertical velocity component, as a 
function of the density variation due to salinity. In 
addition, Gordon and Spaulding (74) have developed a 
three-dimensional hydrodynamics model employing the same 
non-dimensional vertical coordinate used in this study. 

This modal has been designed to interface directly with 
the three-dimensional water quality model. The coupling 
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of these two techniques should provide an effective tool 
for modeling complex pollutant transport in the coastal 
2one. 

Although the information used in the Providence 
River application produced a reasonably good 
verification, it is hoped to have more field information 
than was available for this study. Vertical 
concentration profiles at river mouths, at a number of 
different states of the tide, would permit more accurate 
modalirio of the river source levels. Any information on 
the time variation of outfall discharges would increase 
the value of the modeling effort. The lack of 
information on vertical diffusion coefficients requires 
some arbitrariness. Although field studies attempting 
to define D would be difficult and expensive, they 
would also increase the knowledge to be gained from 
modeling, by revealing the magnitude of vertical 
diffusion in the particular area of interest. As 
mentioned, field studies would also aid in the modeling 
of coTi forms or other indicator bacteria. 

Th- more information that is supplied to the model, 
the more it can produce. The computational methods have 
been shown to be accurate, stable, and very flexible. 
Here is a tool which can help fill many gaps in the 
knowledge of coastal transport problems. 
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